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1.  CONTRACT  PURPOSE 

The  purpose  of  this  contract  (F19628-96-C-0149)  was  to  carry  out  a  program  of  instrument 
construction,  technique  development,  data  acquistion,  and  data  analysis  in  order  to  characterize 
the  auroral  electrojet  and  the  ambient  and  modified  D-region  directly  above  and  near  the  HAARP 
(High  Frequency  Active  Auroral  Research  Program)  facility  in  Alaska.  To  accomplish  this  goal, 
five  VLF  D-region  diagnostic  systems  were  constructed  and  deployed  in  five  High  Schools  near 
the  HAARP  facility  and  along  the  west  coast  of  Alaska.  In  addition,  to  characterize  the  auroral 
electrojet  on  a  continental  scale  and  to  allow  prediction  of  its  location,  four  VLF  D-region 
diagnostic  systems  were  constructed  and  deployed  in  four  High  Schools  on  the  east  coast  of 
the  United  States  and  Canada.  Data  from  the  complete  array  of  D-region  diagnostic  systems 
was  acquired  during  a  number  of  Fall  and  Spring  HAARP  campaigns  during  the  period  of 
performance.  This  data  was  then  used  to  determine  the  state  of  the  D-region  above  HAARP 
during  the  campaigns  and  to  provide  information  on  the  location  of  the  auroral  electrojet.  This 
information  was  essential  to  the  understanding  of  how  the  HAARP  HF  signal  interacted  with 
the  charged  constituents  of  the  D-region  to  produce  ULF/ELFAT.F  waves. 

2.  WORK  PROVIDED 

During  the  period  of  performance  of  this  contract,  Stanford  University  carried  out  the  following 
tasks: 

Task  1 :  Stanford  performed  contract  management,  reporting,  and  technical  overview  activities 
of  contracted  work. 

Task  2:  Stanford  constructed  and  assembled  nine  identical  VLF  D-region  Diagnostics  Systems 
whose  purpose  was  to  allow  the  characterization  of  the  auroral  electrojet  and  the  ambient  and 
modified  D-region  directly  above  and  near  the  HAARP  facility  in  Alaska. 

Task  3:  Stanford  identified  nine  High  Schools  in  which  to  house  the  VLF  D-region  diagnostics 
systems,  and  erected  loop  antennas  at  these  sites  to  detect  three  subionospheric  VLF  signals  from 
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three  VLF  transmitters  that  propagated  above  the  HAARP  facility.  At  each  site  the  system  was 
linked  to  the  internet  with  the  help  of  the  local  High  School  science  teacher. 

Task  4:  Stanford  supplied  computer  code  to  digitally  measure  at  each  High  School  the  phase 
and  amplitude  of  three  subionospheric  VLF  signals  from  three  VLF  transmitters  which  prop¬ 
agated  above  the  HAARP  facility.  These  data  files  were  sent  through  the  internet  to  Stanford 
where  the  data  were  analyzed  and  posted  on  a  Stanford  internet  page  which  was  open  to  the 
public  as  well  as  HAARP  scientists  . 

Task  5:  Stanford  acquired  data  with  the  deployed  array  of  VLFD-region  diagnostic  systems  in 
support  of  each  Spring  and  Fall  HAARP  campaign  during  the  period  of  performance.  This  data 
was  analyzed  and  the  scientific  results  were  presented  at  the  annual  RF  Ionospheric  Interactions 
Workshop  m  Santa  Fe,  New  Mexico.  Experiment  results  were  also  published  in  the  scientific 
literature,  and  formed  one  of  the  bases  of  a  Ph.D.  dissertation  at  Stanford  University. 

3.  INSTRUMENTATION 

Each  of  the  nine  VLF  D-region  Diagnostics  Systems  constructed  under  this  contract  operates  in 
the  same  manner.  At  each  site  the  VLF  signals  from  the  loop  antenna  is  fed  into  a  preamplifier 
box  and  sent  to  a  line  receiver  where  the  signals  are  broadband  filtered  to  remove  high  frqeuency 
noise  and  clipped  to  prevent  saturation  from  nearby  lightning  sources.  The  signal  is  sampled  by 
an  A/D  network  under  computer  control  at  a  rate  of  100  kHz,  and  the  data  is  buffered  to  a  hard 
disk.  Accurate  time  is  provided  through  a  GPS  receiver.  The  buffered  data  are  then  formatted 
and  sent  to  a  Stanford  University  server  by  FTP  over  the  internet. 

Data  from  the  hard  disk  are  used  in  one  of  two  ways.  For  determining  the  phase  and  amplitude 
of  subionospheric  VLF  signals,  the  data  are  digitally  filtered  to  provide  a  pass  band  of  approxi¬ 
mately  100  Hz  centered  on  the  VLF  carrier  of  interest.  This  low  bit  rate  data  is  then  analyzed  to 
obtain  signal  phase  and  amplitude.  To  determine  the  character  and  distribution  of  VLF  signals 
over  the  entire  50  kHz  passband,  the  data  is  Fourier  analyzed  and  presented  as  a  spectrogram. 
This  high  bit  rate  data  is  important  in  identifying  local  sources  of  VLF  noise  that  may  affect  the 
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narrow  band  data,  but  requires  much  more  memory.  Thus  it  is  generally  used  only  in  special 
circumstances. 

4.  RESULTS  OF  HAARP  CAMPAIGNS 

Data  acquired  with  the  VLF  D-region  Diagnostics  Systems  during  the  period  of  performance 
have  been  analyzed  and  the  results  reported  at  the  annual  RP  Ionospheric  Interactions  Workshops 
in  Santa  Fe,  New  Mexico.  Results  have  also  been  published  in  the  scientific  literature  and  used 
in  a  Ph.D.  dissertation  at  Stanford  University,  as  documented  below.  This  Ph.D.  dissertation  is 
reproduced  in  the  Appendix. 

Inan,  U.  S.,  N.  G.  Lehtinen,  S.  J.  Lev-Tov,  M.  P.  Johnson,  T.  F.  Bell,  and  K.  Hurley,  Ionization 
of  the  lower  ionosphere  by  7  -rays  from  a  magnetar:  detection  of  a  low  energy  (3-10  keV) 
component,  Geophys.  Res.  Lett.,  26,  3357,  1999. 

Demirkol,  M.  K.,  U.  S.  Inan,  T.  F.  Bell,  S.  G.  Kanekal,  and  D.  C.  Wilkinson,  Ionospheric 
effects  of  relativistic  electron  enhancement  events,  Geophys.  Res.  Lett.  26,  3557,  1999. 

Demirkol,  M.  K.,  VLF  remote  sensing  of  the  ambient  and  modified  lower  ionosphere,  Ph.D. 
Dissertation,  Stanford  University,  November,  1999. 
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5.  APPENDIX 

The  methodology  involved  in  determining  the  modified  D-region  temperature  and  collision 
frequency  profiles  utilizing  the  data  collected  at  the  three  HAARP  VLF  D-region  Diagnostics 
sites  is  described  in  the  Ph.D.  thesis  of  Dr.  M.  K.  Demirkol.  This  thesis  is  presented  on  the 
following  pages. 
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Abstract 


Electron  density  and  temperature  changes  in  the  D  region  are  sensitively  manifested 
as  changes  in  the  amplitude  and  phase  of  subionospheric  Very  Low  Frequency  (VLF) 
signals  propagating  beneath  the  perturbed  region.  Both  localized  and  large  scale 
disturbances  (either  in  electron  density  or  temperature)  in  the  D  region  cause  signifi¬ 
cant  scattering  of  VLF  waves  propagating  in  the  earth-ionosphere  waveguide,  leading 
to  measurable  changes  in  the  amplitude  and  phase  of  the  VLF  waves.  Large  scale 
auroral  disturbances,  associated  with  intensification  of  the  auroral  electrojet,  as  well 
as  ionospheric  disturbances  produced  during  relativistic  electron  enhancements,  cause 
characteristic  changes  over  relatively  long  time  scales  that  allow  the  assessment  of  the 
‘ambient’  ionosphere.  Localized  ionospheric  disturbances  are  also  produced  by  pow¬ 
erful  HF  transmitting  facilities  such  as  the  High  Power  Auroral  Stimulation  (HIPAS) 
facility,  the  High  frequency  Active  Auroral  Research  Program  (HAARP),  and  also  by 
lightning  discharges.  Amplitude  and  phase  changes  of  VLF  waveguide  signals  scat¬ 
tered  from  such  artificially  heated  ionospheric  patches  are  known  to  be  detectable. 
In  this  study,  we  describe  a  new  inversion  algorithm  to  determine  altitude  profiles  of 
electron  density  and  collision  frequency  within  such  a  localized  disturbance  by  using 
the  measured  amplitude  and  phase  of  three  different  VLF  signals  at  three  separate 
receiving  sites.  For  this  purpose  a  new  optimization  algorithm  is  developed  which 
is  primarily  based  on  the  recursive  usage  of  the  three  dimensional  version  of  the 
Long  Wave  Propagation  Capability  (LWPC)  code  used  to  model  the  subionospheric 
propagation  and  scattering  of  VLF  signals  in  the  earth-ionosphere  waveguide  in  the 
presence  of  ionospheric  disturbances. 
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Chapter  1 
Introduction 


This  dissertation  describes  the  development  of  an  inversion  process  which  uses  ampli¬ 
tude  and  phase  changes  of  coherent  subionospheric  electromagnetic  signals  radiated 
by  very-low-frequency  transmitters  to  determine  the  electron  collision  frequency  pro¬ 
file  in  a  localized  region  of  the  Earth’s  lower  ionosphere  which  is  artificially  modified 
by  powerful  High  Frequency  radiowave  heating.  The  Very-Low-Frequency  (VLF)  sig¬ 
nals  constitute  the  portion  of  the  spectrum  between  3  kHz  and  30  kHz,  although  in 
this  work  we  are  mostly  interested  in  the  range  20  kHz  and  30  kHz,  i.e. ,  the  upper  end 
of  the  VLF  spectrum.  The  ionosphere  consists  of  a  number  of  ionized  regions  above 
the  Earth’s  surface,  and  plays  a  very  important  part  in  global  communication  and 
navigation  systems  [Kelley,  1989].  The  so-called  D  region,  lying  between  60  to  90  km 
altitude,  is  a  good  reflector  for  VLF  waves,  the  approximate  height  of  reflection  for 
20-30  kHz  waves  varying  from  ~65  km  during  the  day  to  ~85  km  during  the  night. 
This  lowest  portion  of  the  ionosphere  is  the  partially  ionized  atmospheric  region  that 
is  generally  a  difficult  region  to  measure,  being  too  high  for  aircraft  and  for  balloons, 
too  low  for  orbiting  satellites,  and  having  too  few  electrons  for  effective  measurements 
by  high  and  very  high  frequency  direct  line-of-sight  radar  systems.  Ground-based  ra¬ 
dio  remote  sensing  at  very  low  frequencies  is  the  most  effective  method  for  studying 
this  region  [Sechrist,  1974]. 

The  fact  that  VLF  waves  are  effectively  reflected  by  the  D  region  can  be  exploited 
to  measure  the  electrical  conductivity  of  this  region.  The  D  region  conductivity 
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depends  upon  both  the  number  density  and  temperature  of  the  ambient  electrons. 
Thus,  the  amplitude  and  phase  of  VLF  waves  propagating  via  reflection  bet-ween  the 
ionosphere  and  the  conducting  Earth  are  sensitively  dependent  on  changes  in  the 
volume  density  and  temperature  of  the  ambient  electrons.  Effects  of  the  electron 
temperature  on  conductivity  are  typically  expressed  in  terms  of  the  frequency  of  the 
collisions  of  the  electrons  with  neutral  atmospheric  constituents.  The  specific  problem 
investigated  in  this  dissertation  is  the  inversion  of  the  amplitude  and  phase  changes  in 
the  VLF  signals  in  order  to  determine  the  electron  collision  frequency  profile  within 
an  ionospheric  region  illuminated  by  a  powerful  high  frequency  radio  transmitter.  In 
this  chapter  we  provide  a  brief  background  on  the  ionosphere,  the  earth-ionosphere 
waveguide,  very  low  frequency  propagation,  and  the  so-called  Luxembourg  effect. 
We  also  describe  and  formulate  the  scientific  problem  addressed  herein  and  list  the 
contributions  of  this  dissertation. 

1.1  The  Ionosphere 

The  ionosphere  is  a  region  of  the  atmosphere  extending  from  an  altitude  of  about  60 
to  more  than  1000  km  in  which  there  exists  a  significant  number  of  free  electrons  and 
ions.  The  presence  of  these  freely  mobile  carriers  of  electric  current  makes  this  region 
an  electrical  conductor,  which  reflects  waves  over  a  broad  range  of  frequencies. 

The  existence  of  an  ionized  layer  in  the  upper  atmosphere  was  appreciated  at 
the  beginning  of  the  20**’  century.  Marconi  demonstrated  that  radio  waves  could 
propagate  across  large  distances  beyond  the  horizon,  as  if  they  were  guided  between  a 
conducting  layer  and  the  ground.  Even  as  early  as  1878,  Stewart  [Stewart,  1878],  from 
the  studies  of  diurnal  oscillations  in  the  Earth’s  magnetic  field,  indicated  that  these 
variations  might  be  due  to  tidal  oscillations  in  the  upper  atmosphere.  However  only 
in  1924  was  the  existence  of  an  ionized  conducting  layer  fully  proven.  Appleton  and 
Barnett  [Appleton  and  Barnett,  1925]  in  England  and  Breit  and  Thve  [Breit  and  Tuve, 
1926]  in  the  USA,  on  both  sides  of  Atlantic,  conducted  experiments  demonstrating 
the  existence  of  the  ionosphere. 
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The  existence  of  the  ionosphere  is  a  consequence  of  two  commonly  known  proper¬ 
ties  of  the  atmosphere.  First,  the  density  of  the  atmosphere  decreases  with  altitude 
above  the  Earth’s  surface,  and  second,  the  atmosphere  absorbs  much  of  the  ultravio¬ 
let  (UV)  and  higher  energy  radiation  incident  upon  it  from  the  Sun  and  other  cosmic 
sources.  Incoming  high  energy  radiation  ionizes  the  atmospheric  constituents.  At 
heights  of  the  order  of  1000  km  or  more  the  medium  is  almost  fully  ionized,  but  so 
tenuous  that  the  ion  and  the  electron  concentrations  are  relatively  small.  At  lower 
levels  there  is  more  air  to  be  ionized  so  that  the  radiation  is  more  strongly  absorbed 
and  the  electron  and  ion  concentrations  are  higher.  The  ions  and  free  electrons  cre¬ 
ated  this  way  are  absorbed  by  chemical  interactions  (recombination  and  attachment 
processes)  with  other  ions  and  atoms,  leading  to  a  dynamic  equilibrium  between  the 
production  of  the  electrons  and  their  absorption,  resulting  in  an  average  finite  free 
electron  and  ion  density.  At  still  lower  altitudes  the  air  density  is  so  high  that  the 
absorption  process  dominates  and  resulting  in  much  lower  free  electron  densities. 

The  ionosphere  can  be  viewed  as  a  variable  shell  of  plasma  surrounding  the  Earth. 
Above  1000  km  the  atmosphere  is  a  fully  ionized  plasma  and  is  therefore  a  relatively 
good  electrical  conductor.  Furthermore,  the  density  of  air  is  so  low  that  free  electrons 
and  ions  do  not  undergo  collisions  with  other  ions  or  molecules,  and  can  freely  move 
under  the  influence  of  electromagnetic  fields.  Under  these  conditions  the  Earth’s 
magnetic  field  has  a  strong  effect  on  the  structure  of  the  medium,  and  the  region  is 
thus  known  as  magnetosphere.  It  extends  out  to  a  distance  of  14i?E  on  the  sunlit  side 
of  the  Earth,  where  R^  ~  6400  km  is  the  radius  of  the  Earth. 

Historically  the  ionosphere  has  been  divided  into  layers,  the  earliest  detected  being 
the  E  layer,  so  named  due  to  reflection  of  electric  fields.  The  D  and  F  layers  axe 
below  and  above  the  E  layer  respectively.  The  D  region,  between  60  and  90  km  is 
highly  variable  with  much  lower  electron  density.  However,  the  electron  density  in 
all  of  these  regions  widely  varies  with  time  of  day,  season,  solar  cycle  and  the  level 
of  magnetospheric  and/or  solar  wind  activity  [Davies  1965;  p.  265].  The  focus  of 
our  attention  in  this  work  is  the  D  region  of  the  ionosphere,  lying  between  ~60  and 
90  km.  Typical  profiles  of  electron  density  and  identification  of  the  various  altitude 
ranges  are  given  in  Figure  1.1. 
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Figure  1.1:  Typical  midlatitude  ionospheric  electron  density  profile  for  sunspot  maxi¬ 
mum  and  minimum  conditions  during  day-time  and  night-time.  The  different  altitude 
regions  in  the  ionosphere  are  labelled  with  D,  E,  Fi,  and  Fz  [Richmond,  1987], 
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1.2  The  Earth-Ionosphere  Waveguide 

One  of  the  properties  of  electromagnetic  waves  is  that  they  reflect  from  conducting 
boundaries  and  can  thus  be  channelled  or  guided  between  such  boundaries.  The 
surface  of  the  Earth  is  a  moderately  good  conductor  of  electricity,  and  can  reflect 
radio  signals  in  the  lower  frequency  signal  ranges,  including  VLF  waves  [Davies, 
1956:  p.  45].  The  ionosphere,  on  the  other  hand,  is  a  complicated,  highly  lossy 
(conducting),  and  anisotropic  medium  and  the  D  region  is  a  good  reflector  of  VLF 
■waves,  as  mentioned  above.  The  Earth  and  the  ionosphere  thus  constitute  the  two 
boundaries  of  a  waveguide  vdthin  which  VLF  radio  waves  can  propagate,  commonly 
referred  to  as  the  earth-ionosphere  waveguide. 


Figure  1.2:  VLF  signal  propagation  can  be  thought  of  as  repeated  reflection  back  and 
forth  in  the  earth-ionosphere  waveguide. 

The  propagation  of  VLF  electromagnetic  signals  in  the  earth-ionosphere  waveg¬ 
uide  can  be  thought  of  in  terms  of  signals  reflecting  back  and  forth  between  the  Earth 
and  ionosphere  as  shown  in  Figure  1.2,  although  analysis  of  this  propagation  can  be 
more  effectively  described  in  terms  of  waveguide  modes  as  described  in  Section  1.3. 
An  electromagnetic  wave  totally  reflects  from  a  medium  with  \^.rying  dielectric  prop¬ 
erties  at  the  point  at  which  the  refractive  index  is  zero.  For  an  isotropic  lossless 
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plasma,  this  condition  is  realized  when  uj  Ci:  ujp,  where  Up  is  the  plasma  frequency  of 
the  medium  given  by  LOp  =  AirNeq^f  (come) ,  with  N^  being  the  number  density  of  elec¬ 
trons,  Qe  being  the  charge  of  an  electron,  cq  being  the  permittivity  of  free  space,  and 
TUe  being  the  mass  of  an  electron.  However,  in  the  VLF  range,  the  absorption,  reflec¬ 
tion  and  transmission  of  a  radio  signal  incident  on  the  ionosphere  depends  on  factors 
such  as  the  wave  frequency,  the  angle  of  incidence,  the  altitude  profiles  of  electron 
and  ion  concentrations,  the  altitude  profile  of  electron  temperature  and  thus  the  rate 
at  which  these  constituents  collide  with  neutrals  (i.e.,  the  collision  frequency),  and 
the  intensity  and  local  orientation  of  the  Earth’s  magnetic  field.  Due  to  the  influence 
of  the  Earth’s  magnetic  field,  and  the  presence  of  absorption  due  to  collisions,  the 
refractive  index  does  not  reach  zero  at  any  D  region  altitude,  so  total  reflection  never 
occurs,  but  substantial  partial  reflection  is  possible.  The  region  of  the  ionosphere 
where  such  partial  reflection  occurs  can  be  thought  of  as  a  region  where  the  refrac¬ 
tive  index  is  changes  very  rapidly  over  distances  comparable  to  a  wavelength.  WTien 
this  occurs,  the  region  acts  like  a  sharp  boundary  between  two  media  and  reflection 
occurs.  Quantitatively,  the  effective  reflection  height  under  nighttime  conditions  for 
VLF  waves  is  determined  by  the  condition  [Wait  and  Spies,  1964]: 

uP- 

w  ~  ~  2.5  X  10^  s~^  (1.1) 

where  ujp  is  the  plasma  frequency  of  the  medium  given  by  Up  =  ^■nNeql/{eQme),  with 
Ne  being  the  number  density  of  electrons,  being  the  charge  of  an  electron,  eo  being 
the  permittivity  of  free  space,  and  me  being  the  mass  of  an  electron;  uj  is  the  wave 
frequency;  and  u  is  effective  collision  frequency  of  electrons  with  heavy  particles. 

The  reflecting  properties  of  the  ionosphere  are  thus  dependent  on  the  number 
density  of  electrons,  Ne,  which  varies  with  height  above  the  Earth’s  surface  (Figure 
1.3)  and  which  is  particularly  variable  in  the  nighttime  D  region.  A  tvpical  D  region 
electron  density  profile  is  shown  in  Figure  1.3-a  and  in  the  critical  altitude  range  of 
60  to  90  km  can  be  approximately  described  by  a  two-parameter  exponential  profile. 

Ne{h)  =  1.43  X  10^^  exp  [{0  -  0.15)  {h  -  h')]  (1.2) 
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with  Ne  in  electrons/cm^,  h'  being  the  effective  reflection  height  in  km,  and  p  in 
determining  the  sharpness  (slope)  of  the  profile.  Typical  values  for  h  and  P  for 
nighttime  conditions  are  h'  ^  85  km  and  P  c:i  0.5  km-\  respectively  [Wait  and  Spies, 
1964].  The  overall  electron  density  profile  shown  in  Figure  1.3-a  has  been  used  in  a 
range  of  recent  work  involving  modeling  of  VLF  scattering  from  localized  ionospheric 
disturbances  [Inan  et  ai,  1992;  Rodriguez  et  at,  1994;  Lev-Tov  et  al,  1995]. 

Absorption  of  radio  waves  in  the  ionosphere  occurs  primarily  due  to  collisions 
between  electrons  and  the  neutral  constituents  of  the  atmosphere.  Electrons  move 
in  response  to  the  electromagnetic  fields  of  the  propagating  waves  and  collide  with 
the  neutrals  and  ambient  ions.  These  collisions  transfer  the  propagating  wave  energy 
to  thermal  energy,  and  thus  reduce  the  wave  intensity.  The  lower  regions  of  the 
ionosphere  have  the  greatest  density  of  neutrals  and  thus  have  the  highest  rate  of 
occurrence  of  these  collisions  [Ratcliffe,  1972].  A  typical  collision  frequency  profile  as 

a  function  of  altitude  is  shown  in  Figure  1.3-b. 

The  average  number  of  collisions  v  which  an  electron  makes  per  unit  time  with 
the  air  molecules  depends  upon  the  number  density  of  the  molecules,  and  therefore 
is  determined  largely  by  the  density  and  composition  of  the  air  and  also  to  a  lesser 
degree  by  the  velocity  distribution  of  the  electrons  [Budden,  1961].  Various  collision 
firequency  altitude  profiles  exist  for  electron-neutral  particle  collisions  iy(.{z)  [Morfitt 
and  Shellman,  1976].  An  exponential  profile  as  given  in  equation  (1.3)  is  a  typical 
example: 

i/e(z)  =  1.816  X  (1-3) 

where  u  is  in  s-^and  z  is  in  km.  This  profile  is  adopted  as  the  ambient  (unperturbed) 
collision  frequency  profile  in  this  work. 


1.3  VLF  Propagation  Models 

Understanding  and  accurate  modeling  of  the  propagation  of  electromagnetic  waves  in 
the  earth-ionosphere  waveguide  is  essential  for  the  purpose  of  this  dissertation.  Any 
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Figure  1.3:  (a)  Typical  nighttime  electron  density  as  a  function  of  altitude  in  the 
D  region,  (b)  A  typical  model  for  the  effective  electron  collision  frequency  versus 
altitude  in  the  lower  ionosphere. 
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model  of  electromagnetic  wave  propagation  must  be  based  on  a  simple  set  of  formulas 
knows  as  Maxwell’s  equations,  as  given  below  in  differential  form: 


-  dB 

(1.4) 

-  -  dD 

V  xH  =  J+-^ 
ot 

(1.5) 

Ci. 

II 

t> 

(1.6) 

V-5  =  0 

(1.7) 

where 

E  is  the  total  electric  field, 

H  is  the  total  magnetic  field, 

D  =  eE, 

B  =  ijlH, 

e  is  the  permittivity  of  the  medium, 
fjL  is  the  permeability  of  the  medium, 
p  is  the  volume  charge  density,  and 
J  is  the  current  density. 

The  electric  and  magnetic  fields  everywhere  within  the  waveguide  are  obtained  by 
solving  these  equations  subject  to  boundary  conditions.  The  radio  receivers  used  to 
acquire  the  data  utilized  in  this  dissertation  are  designed  to  measure  the  electromag¬ 
netic  fields  at  ground  level. 

There  are  mainly  two  methods  for  studying  the  propagation  of  VLF  signals  in 
the  earth-ionosphere  waveguide.  The  ray  optics  approach  [Bom  and  Wolf,  1965]  is 
generally  easier  to  visualize  and  analyze.  This  method  simply  traces  the  different 
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discrete  ray  paths  that  complete  an  integral  number  of  reflections  between  the  trans¬ 
mitter  and  the  receiver  (Figure  1.4).  The  number  of  complete  reflections  of  the  wave 
between  the  Earth  and  ionosphere  before  it  reaches  the  receiver  is  called  the  order 
of  the  ray  path.  As  the  order  of  the  ray  path  increases  the  angle  of  incidence  on  the 
ionosphere  also  increases.  Since  the  absorption  and  transmission  rates  are  higher  for 
higher  angles  of  incidence,  the  total  effect  of  these  signals  at  the  receiver  becomes  less 
and  less  significant  with  increasing  angle  of  incidence. 

In  ray  theory  few  ray  paths  need  to  be  calculated  for  very  short  (<500)  km  prop¬ 
agation  paths.  However  when  the  propagation  path  is  longer,  there  are  more  ray 
paths  with  low  attenuation  between  the  transmitter  and  receiver  and  calculations  be¬ 
come  tedious.  In  addition,  for  each  separation  distance  between  the  transmitter  and 
receiver  the  ray  path  geometry  changes  and  all  of  the  calculations  must  be  repeated. 

Electromagnetic  waves  propagating  in  the  earth-ionosphere  waveguide  can  also  be 
represented  in  terms  of  a  sum  of  an  infinite  series  of  discrete  waveguide  modes.  The 
order  of  a  waveguide  mode  simply  corresponds  to  the  number  of  maxima  and  minima 
in  the  transverse  vertical  field  pattern.  As  the  waveguide  mode  order  increases,  the 
spatial  attenuation  factor  increases.  Thus,  the  higher  order  modes  do  not  propagate 
as  far  in  the  waveguide  as  the  lower  order  modes. 

Each  waveguide  mode  can  be  represented  in  terms  of  a  superposition  of  two  uni¬ 
form  plane  waves  propagating  by  reflecting  back  and  forth  between  the  waveguide 
walls,  being  incident  on  the  boundaries  at  a  discrete  angle  from  the  vertical  [Inan 
and  Inan,  1999;  p.  269).  The  discrete  set  of  complex  angles  9n  for  which  this  occurs 
are  called  the  eigenangles,  each  of  which  corresponding  to  a  single  waveguide  mode. 

Accurate  modeling  of  propagation  in  the  earth-ionosphere  waveguide  is  not  com¬ 
plete  if  the  influence  of  the  Earth  curvature  is  not  accounted  for.  If  this  approach  is 
used  to  yield  an  electric  field  E  and  a  magnetic  field  H  whose  components  are  each 
functions  of  'F  {r,9),  an  infinite  sum  of  rather  complicated  spherical  wave  functions 
[Wait,  1970;  p.  157]  results  ; 


{r,G)  =  '^DqZq{kr)Pq{-co^6)  (1.8) 

n=l 
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Figure  1.4:  (a)  Ray  Paths.  Possible  discrete  raj^  paths  from  a  transmitter  to  a 
receiver.  There  are  an  integral  number  of  reflections.  Only  the  first  three  ray  paths  are 
indicated.  For  the  n^'"  order  ray  path  there  are  2”~^ reflection  points  in  tire  ionosphere, 
(b)  Waveguide  modes.  The  vertical  \^iation  of  the  field  amplitude  of  different 
modes  are  shown.  The  labels  Fq  and  Fi  are  referred  to  later  in  connection  with 
equation  (1.12).  Also  shown  are  ra^aths  for  component  uniform  plane  waves  (solid 
and  dashed)  the  superposition  of  which  can  be  used  to  represent  the  waveguide  modes. 
Each  mode  corresponds  to  an  incidence  angle  dn  for  the  component  waves. 
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where 

r,  0  are  spherical  coordinates, 
n  is  the  mode  number, 

5  is  a  complex  number, 
are  the  mode  coefficients, 

are  Hankel  functions  of  the  first  and  second  kind  of  order  i/, 
k  =  uj/c, 
u  =  27r/, 

/  is  the  wave  frequency,  and 
Pt,  are  Legendre  functions. 

It  was  shown  1970;  p.  136]  that  under  undisturbed  ionospheric  conditions, 

with  the  space  between  the  Earth  and  the  ionosphere  taken  to  constitute  a  spherically 
concentric  waveguide,  and  with  homogeneous  conditions  along  the  entire  path,  the 
total  electric  field  Etotai  s-t  a  great  circle  path  (GCP)  distance  d  from  the  transmitter 
can  be  written  as: 


1  °° 

“  VIsin  WflE)l  (1-9) 

where 

n  is  the  mode  number, 
i  = 

6n  is  the  complex  valued  angle  of  incidence  of  each  mode  on  the  ionosphere,  known 
as  the  eigenangle, 

S„  =  sin(0„),  is  the  sine  of  the  complex  eigenangle  On, 
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Re  is  the  radius  of  Earth, 

=  KTnRn, 

K  =  a^/Pj,  where 
a  is  a  constant  coefficient, 

P  is  the  radiated  power, 

/  is  the  signal  frequency,  and  where 


(1.10) 
(1.11) 

j  =  1,2,3,  represent  the  three  dimensions,  x,y,z  where  the  GCP  lies  along  the  x 
direction, 

tj  and  Tj  are  factors  for  the  orientation  of  the  transmitter  and  the  receiver  antenna 
respectively, 

and  A^„  are  the  initial  excitation  factors  for  each  mode  with  respect  to  the 

JtiTI 

transmitter  and  receiver  antenna  respectively, 

Glni^)  and  C?f„(z)  are  the  height  gain  factors  for  each  mode  evaluated  at  the 
transmitter  and  receiver  antenna  respectively. 

When  the  waveguide  walls  are  not  homogeneous  but  exhibit  slowly  varying  con¬ 
ductivity,  the  angle  of  incidence  is  a  function  of  the  propagation  path  and  the 
sin(0„)  term  in  equation  (1.9)  is  replaced  by  j*  sin  [0„  (x)]  dx. 

Any  quantitative  use  of  equation  (1.9)  requires  that  the  complex  eigenangles  ^ 
are  known.  The  eigenangles  9n  can  be  determined  using  Budden’s  formulation  of  the 


3 

j  =  l 


33 


CHAPTER  1.  INTRODUCTION 


14 


earth-ionosphere  waveguide  problem  [Budden,  1961;  p.  115-123].  In  this  formulation, 
a  modal  equation  is  developed  in  terms  of  reflection  coefficients  at  the  ionosphere 
and  Earth,  designated  respectively  as  and  where  is  the  angle  of 

incidence  of  the  waves  on  the  ionosphere  as  above.  The  ionosphere  is  an  anisotropic 
medium  due  to  the  existence  of  Earth’s  magnetic  field,  thus  the  reflection  coefficient 
is  actually  a  coefficient  matrix.  Since  the  medium  is  lossy,  the  reflection  matrix  is 
complex.  In  order  to  have  the  upgoing  wave  and  the  twice  reflected  wave  interfere 
constructively,  the  twice  reflected  wave  must  have  a  total  phase  change  of  27rm.  Thus 
if  we  let  Fq  denote  the  upgoing  wave  and  the  F\  denote  the  twice  reflected  wave  (See 
Figure  1.4),  the  coherence  or  self-consistent  mode  condition  is  given  in  (1.12). 

The  eigenangles  are  determined  from  this  equation  by  setting  F\  =  Foj . 

Fl|  R9  j2) 

At  the  simplest  level,  we  can  approximate  the  earth-ionosphere  waveguide  by 
treating  the  ground  as  a  perfect  electric  conductor  and  the  ionosphere  as  a  perfect 
magnetic  conductor  [Davies,  1965;  p.  423].  In  this  case  F'  =  -1,  F®  =  1,  and  the 
resultant  eigenangles  obtained  from  equation  (1.12)  are  shown  in  Figure  1.5. 

1.4  Ionospheric  Heating 

The  basic  thrust  of  this  dissertation  is  to  develop  and  apply  a  new  method  of  diag¬ 
nostics  of  the  lower  ionosphere  under  conditions  when  it  is  modified  by  a  powerful  HF 
heater.  The  method  utilized  relies  upon  the  measurement  of  the  modulation  pattern 
of  the  HF  heater  as  it  is  imposed  on  a  VLF  radio  beacon  signal.  The  fundamental 
physical  means  of  coupling  of  the  HF  and  VLF  signals  is  the  nonlinear  response  of 
the  ionosphere  to  the  powerful  HF  heater  wave.  One  of  the  first  recognized  manifes¬ 
tations  of  ionospheric  heating  by  radio  waves  was  the  so-called  ’Luxembourg  effect’, 
involving  such  nonlinear  coupling  between  two  radio  signals. 
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Real(ej^)[deg] 


Figure  1.5:  When  the  Earth  surface  is  treated  as  a  perfect  electric  conductor  and  the 
ionosphere  is  treated  as  a  perfect  magnetic  conductor,  we  find  R^  =  1  and  R^  =  —  1. 
Thus  the  resonance  mode  condition  gives  us  the  eigenangles  in  the  figure  for  h  =  85 
km  and  k  =  w/c—  (27r)(25  kHz)/(3  x  10®m/s) =5.2359  x  10“^m"^  The  numbers  refer 
to  the  mode  numbers  n. 
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1.4.1  The  Luxembourg  Effect 

In  1933,  it  was  discovered  that  the  AM  radio  transmissions  received  in  Eindhoven, 
Netherlands  from  a  station  in  Beromunster,  Switzerland,  frequently  had  superimposed 
upon  them  the  modulation  of  transmissions  at  a  different  radio  frequency  broad¬ 
cast  from  Luxemburg  [Tellegen,  1933].  This  event,  called  the  Luxembourg  effect  or 
ionospheric  cross-modulation  or  wave  interaction,  was  interpreted  to  be  caused  by 
a  nonlinearity  in  the  ionosphere.  The  geometry  involved  in  the  ionospheric  cross 
modulation  is  illustrated  in  Figure  1.6. 

Cross-modulating  region 


Figure  1.6.  Geometry  and  nomenclature  involved  in  ionospheric  cross-modulation. 

A  receiver  is  tuned  to  receive  a  signal  called  the  wanted  wave  ( shown  as  initially 
unmodulated)  from  a  wanted  transmitter.  The  received  signal  shows  modulation, 
called  transfer  modulation,  imposed  on  the  wanted  wave  in  the  cross-modulating  region 
under  the  action  of  the  disturbing  wave.  ^Vhen  a  powerful  disturbing  wave  is  absorbed 
in  the  cross-modulation  region,  this  absorbed  energy  goes  partially  into  an  increase 
in  the  mean  energy  of  the  electrons,  increasing  the  electron  collision  frequency.  The 
modulation  of  the  disturbing  wave  therefore  produces  corresponding  variations  in  the 
collision  frequency  of  the  region.  Consequently,  the  absorption  of  the  transmitted  and 
wanted  wave  which  propagates  through  this  region  exhibits  the  same  time  variation; 
i.e.,  the  modulation  is  “transferred”  from  the  disturbing  wave  to  the  wanted  wave 
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[Bailey  and  Martyn,  1934]. 

The  initial  studies  of  ionospheric  modification  in  the  context  of  Luxembourg  ef¬ 
fect  have  evolved  into  a  new  field  of  high  frequency  (HF)  ionospheric  modification 
studies  conducted  with  dedicated  experimental  heating  transmitters.  Powerful  HF 
heating  facilities  have  been  constructed  in  the  United  States,  Norway,  and  Russia  for 
controlled  ionospheric  experiments  [Wang  and  Brandt,  1990;  and  references  therein]. 
Although  these  ionospheric  heating  experiments  have  a  wide  range  of  objectives,  one 
of  the  documented  effects  of  heating  of  the  ionosphere  by  HF  waves  is  the  resultant 
enhancement  of  increasing  the  electron  collision  frequency  in  the  D  region  by  an  or¬ 
der  of  magnitude  over  areas  greater  than  300  km^  [Stuhhe  and  Kopka,  1979].  Such 
large  changes  in  the  D  region  in  turn  affect  the  propagation  of  VLF  waves  in  earth- 
ionosphere  waveguide  as  has  been  first  observed  by  Barr  et  al.  [1985]  and  as  further 
discussed  in  the  next  subsection. 

1.4.2  HF/VLF  Coupling 

The  feasibility  of  measuring  the  amplitude  and  phase  of  the  VLF  waveguide  signals 
scattered  by  artificially  heated  portions  of  the  ionosphere  has  been  demonstrated  by 
different  experiments  [Barr  et  al,  1985;  Dowden  et  al,  1991].  In  this  section  we  briefly 
describe  an  experiment  which  documented  VLF  signatures  of  ionospheric  heating  by 
the  High  Power  Auroral  Stimulation  (HIPAS)  HF  facility  (located  near  Fairbanks, 
Alaska)  in  1992  [Bell  et  al,  1995]. 

In  the  Fall  1992  and  Spring  1993,  an  ionospheric  modification  campaign  was  car¬ 
ried  out  to  determine  the  characteristics  of  the  lower  ionosphere  during  the  operations 
of  the  HIPAS  HF  heater.  The  VLF  data  described  in  this  section  were  acquired  at 
Fort  Yukon  (FY),  Alaska.  The  HIPAS  facility  operated  with  a  total  radiated  power  of 
800  kW  and  an  effective  radiated  power  of  55  MW  (ERP)  at  2.85  MHz.  The  receiver 
at  FY  continuously  monitored  the  signals  from  VLF  transmitters  NPM  in  Hawaii, 
NLK  in  Washington  and  Silver  Creek  (SC)  in  Nebraska.  The  configuration  of  the 
great  circle  paths  between  the  transmitters  and  the  receiver  is  shown  in  Figure  1.7. 
Table  1.1  shows  latitude  and  longitude  information  for  the  VLF  receiver  site  FY 
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Latitude 

Longitude 

FY  (VLF  Receiver) 

66.6°N 

145.2°W 

HIPAS  (Heating  Facility) 

64.9°N 

146.8°W 

Table  1.1;  Locations  of  VLF  Receiver  FY  and  heating  facility  HIPAS  in  Alaska. 


Transmitter  Name 

Latitude 

Longitude 

frequency 

(kHz) 

Power 

(kW) 

NLK  (Jim  Creek,  Washington) 

48°  12'  N 

121°  55'  W 

24.8 

850 

SC  (Silver  Creek,  Nebraska) 

41°  30'  N 

97°  36'  N 

48.5 

60 

NPM  (Lualualei,  Hawaii) 

2r25'N 

158°  09'  N 

23.4 

300 

Table  1.2:  List  of  VLF  transmitters  used  during  HIPAS  1992  Campaign 


and  the  HIPAS  observatory,  Table  1.2  shows  the  list  of  VLF  transmitters  used  in  this 
experiment. 

The  great  circle  path  of  the  NPM  signal  received  at  FY  passes  within  20  km  of 
HIPAS.  The  distances  between  the  NLK-FY  and  SC-FY  great  circle  paths  and  the 
HIPAS  facility  were  much  larger  and  thus  VLF  signals  from  NLK  and  SC  did  not 
show  any  sign  of  ionospheric  heating.  But  the  amplitude  and  phase  of  the  signal 
from  the  NPM  transmitter  observed  in  Fort  Yukon  exhibited  a  clear  change  with  the 
same  ON/OFF  modulation  pattern  as  that  of  the  HIPAS  HF  transmissions.  For  28 
minutes,  the  HIPAS  vertical  beam  (20°  full  beam  width  at  half  maximum  power) 
was  switched  ON  and  OFF  at  2  Hz.  Figure  1.8  shows  the  results  of  this  experiment. 
The  upper  panel  shows  the  NPM  amplitude  during  the  28  minutes  the  experiment 
was  conducted.  The  middle  left  panel  is  the  superposed  epoch  analysis  which  is 
obtained  simply  by  dividing  the  data  in  the  upper  panel  into  500  ms  serial  sections 
and  subsequently  summing  and  averaging  to  obtain  the  final  single  500  ms  result. 
The  first  100  ms  segment  of  each  500  ms  period  corresponds  to  the  ON  period  of 
the  HF  heater.  When  the  heater  is  OFF,  the  total  signal  received  at  FY  is  simply 
the  unperturbed  NPM  signal  which  propagates  directly  to  FY.  However  when  the 
heater  is  ON,  in  addition  to  the  direct  signal,  there  is  the  signal  scattered  from  the 
heated  ionosphere  over  the  HIPAS  heater,  resulting  in  an  amplitude  increase  of  about 
0.18  dB  in  the  total  received  signal.  In  other  words,  the  100  ms  ON/400  ms  OFF 
modulation  of  the  HF  heater  is  ‘transferred’  to  the  NPM-FY  VLF  signal,  precisely 
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in  the  same  manner  as  the  classical  cross-modulation  problem  depicted  in  Figure  ??. 
The  lower  left  panel  shows  the  spectral  analysis  of  the  raw  data.  Peaks  at  2  Hz  and 
its  harmonics  are  clear  indications  of  the  VLF  modulation  resulting  from  the  pulse 
modulation  of  the  HF  heater  signal.  The  panels  on  the  right  show  the  superposed 
epoch  analysis  of  the  measured  phase  of  NPM  signal  and  the  spectral  analysis  of  the 
raw  phase  data.  The  HF  heating  modulation  is  clearly  evident  on  the  NPM  signal 
phase,  with  the  phase  change  being  ~4.5°. 

1.4.3  VLF  Diagnostics  of  the  Modified  D  region 

The  Stanford  VLF  observations  during  the  1992  HIPAS  campaign  were  conducted 
as  part  of  the  development /testing  of  a  new  diagnostic  technique  whereby  the  phase 
and  amplitude  changes  of  the  VLF  signals  propagating  nearby  a  heated  ionospheric 
region  are  used  to  determine  the  collision  frequency  profile  in  the  heated  region.  This 
technique  was  developed  as  a  diagnostic  tool  for  the  new  HAARP  (HF  Active  Auroral 
Research  Program)  ionospheric  heating  facility,  now  operational  in  Alaska  (Chapter 
5).  In  this  subsection  we  briefly  describe  the  formulation  of  this  new  diagnostic  tool, 
the  various  components  of  which  constitute  the  core  of  this  Ph.D.  dissertation. 

The  measurement  of  any  given  VLF  signal  phase  and  amplitude  changes  produced 
by  a  HF  heater  provides  data  from  which  we  can  characterize  the  electron  collision 
frequency  profile  with  two  free  parameters.  The  measurement  of  the  phase  and  am¬ 
plitude  changes  in  N  distinct  signals  (at  N  different  frequencies)  thus  provides  data 
from  which  we  can  characterize  this  profile  with  2N  free  parameters.  Due  to  practical 
considerations  (e.g.,  availability  of  VLF  transmitters  and  field  deployment  sites),  the 
Stanford  VLF  D  region  diagnostic  technique  is  based  on  the  use  of  3  distinct  VLF 
signals,  as  described  in  the  next  paragraph. 

Amplitude  and  phase  changes  of  three  different  VLF  signals  received  at  three 
different  sites  are  measured  with  high  time  resolution  both  before,  during  and  after 
episodes  of  HF  heating  with  prescribed  modulation  patterns.  The  transferred  modu¬ 
lation  on  the  subionospheric  VLF  signals  are  then  extracted  via  analysis  (similar  to 
that  shown  in  Figure  1.8)  and  are  interpreted  in  terms  of  altitude  profiles  of  D 
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Figure  1.8:  (a)  Raw  NPM  aonplitude  data  recorded  at  FY  while  the  HIPAS  receiver 
w'as  ON  for  100  ms  and  OFF  for  400  ms  (2Hz  peridodicity),  repeated  for  28  min. 
(b)  Superposed  epoch  analysis  (with  T=500  ms)  of  the  raw  data  set  from  (a),  (c) 
Spectral  analysis  of  the  raw  data  set  from  (a),  (d)  Superposed  epoch  analysis  with 
T=500  ms  of  raw  NPM  phase  data  at  FY  during  whi^  the  HIPAS  HF  heater  was 
on  for  100  ms  and  off  for  400  ms  (2  Hz  periodicity)  (e)  Spectral  analysis  of  the  same 
raw  data  [Bell  et  al.,  1995]. 
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region  electron  collision  frequency,  or  equivalently,  electron  temperature  profiles. 

This  VLF  long  path  diagnostic  technique  relies  heavily  on  sophisticated  models 
of  VLF  propagation  and  scattering,  as  indicated  in  Figure  1.9.  The  propagation  of 
the  VLF  signal  in  the  earth-ionosphere  waveguide  along  the  great  circle  path  from 
the  VLF  transmitter  to  the  receiver,  and  its  scattering  from  the  heated  ionospheric 
patch,  are  modeled  with  the  Long  Wave  Propagation  Capability  (LWPC)  code  [Fergu¬ 
son  and  Snyder,  1987]  and  a  three  dimensional  version  of  this  code  recently  developed 
at  Stanford,  both  described  in  detail  in  Chapter  2.  The  use  of  the  VLF  propagation 
and  scattering  codes  for  extraction  of  D  region  electron  collision  frequency  and  tem¬ 
perature  profiles  firom  VLF  amplitude  and  phase  data  from  multiple  sites  as  shown 
in  Figure  1.9  constitutes  the  new  VLF  D  region  diagnostic  tool  which  is  the  subject 
of  this  dissertation. 


Figure  1.9:  Development  of  algorithm  for  extraction  of  D  region  electron  collision 

firequency  profile  from  VLF  amplitude  and  phase  changes  measured  at  three  different 
sites. 


1.5  Contributions 

The  contributions  of  this  dissertations  are  as  follows  : 
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•  Interpretation  of  an  observed  correlation  between  relativistic  electron  precipi¬ 
tation  enhancement  events  and  subionospheric  VLF  signal  amplitude  in  subau- 
roral  regions.  VLF  signal  amplitudes  exhibit  >  10  dB  changes,  exhibiting  the 
same  27  day  cycle  and  2-3  day  rise  and  fall  time  pattern  as  relativistic  electron 
enhancement  events  recorded  by  GOES  7  and  SAMPEX.  This  new  observation 
[Demirkol  et  al,  1999]  indicates  that  the  nighttime  lower  ionospheric  electron 
density  at  subauroral  latitudes  is  detectably  affected  by  27-day  periodicity  in 
solar  rotation.  The  experimental  data  and  its  interpretation  is  presented  in 
Chapter  3. 

•  Development  of  a  method  to  determine  the  ambient  electron  density  profile  over 
a  subauroral  propagation  path  based  on  the  fact  that  D  region  electron  density 
profile  is  strongly  modulated  by  precipitation  enhancement  events.  This  new 
method  is  described  in  Section  3.2. 

•  Development  of  a  new  inversion  algorithm  to  determine  the  altitude  profile  of 
electron  collision  frequency  (or  temperature)  within  localized  disturbances  by 
using  the  measured  amplitude  and  phase  of  three  different  VLF  signals  at  three 
separate  receiving  sites.  Localized  disturbances  are  routinely  produced  by  pow¬ 
erful  HF  transmitting  facilities  such  as  the  High  Power  Auroral  Stimulation 
(HIPAS)  facility  [Bell  et  al,  1995]  and  the  High  frequency  Active  Auroral  Re¬ 
search  Program  (HAARP)  and  also  by  lightning  discharges  [Inan  et  al,  1996]. 
The  new  optimization  algorithm  developed  here  is  primarily  based  on  the  re¬ 
cursive  usage  of  the  three  dimensional  version  [Poulsen  et  al,  1993]  of  the  Long 
Wave  Propagation  Capability  (LWPC)  code  used  to  model  the  subionospheric 
propagation  and  scattering  of  VLF  signals  in  the  earth-ionosphere  waveguide 
in  the  presence  of  ionospheric  disturbances.  The  new  inversion  algorithm  is 
described  and  documented  in  Chapter  4. 

•  Application  of  the  inversion  techniques  to  experimental  data  acquired  during 
the  March  1999  HAARP  campaign.  The  new  inversion  technique  is  used  to 
extract  the  first  measurements  of  the  effective  electron  collision  frequency  in  a 
heated  region,  as  described  in  Chapter  5. 
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Chapter  2 

VLF  Signal  Propagation  Models 


VLF  waves  at  large  distances  (>~1000  km)  from  their  source  can  be  conveniently 
represented  in  terms  of  a  summation  of  multiple  waveguide  modes,  each  propagating 
with  a  different  complex  eigenangle  in  the  waveguide  formed  by  the  conducting  curved 
Earth  and  the  anisotropic,  imperfectly  conducting  curved  ionosphere  with  a  dipping 
magnetic  field.  In  this  chapter,  we  describe  the  tools  we  use  to  quantitatively  analyze 
the  propagation  and  scattering  of  such  signals  and  also  present  an  analysis  of  a  single 
waveguide  mode  case. 

2.1  Long  Wave  Propagation  Capability 

The  Long  Wave  Propagation  Capability  (LWPC)  code,  developed  over  many  years 
by  the  Naval  Oceans  System  Center  [Ferguson  and  Snyder,  1987],  and  now  widely 
used,  is  a  combination  of  FORTRAN  programs  which  enable  the  user  to  model  the 
propagation  of  VLF  signals  in  two  dimensions  along  the  Great  Circle  Path  (GCP)  be¬ 
tween  a  transmitter  and  a  receiver.  Following  is  a  brief  description  of  these  programs 
together  with  their  input  and  output  specifications. 

Given  the  required  parameters,  the  electric  field  along  a  propagation  path  can 
be  calculated  using  the  straightforward  two  dimensional  method  given  by  equation 
(1.9).  In  this  case  the  dimensions  are  x  and  z,  where  x  is  along  the  great  circle 
path  between  the  transmitter  and  the  receiver,  and  z  is  the  altitude,  and  all  of  the 
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properties  of  the  ionosphere  and  ground  are  assumed  to  be  constant  in  the  y  direction. 
Figure  2.1  shows  the  geometry  for  the  2D  \XF  propagation  problem. 


Figure  2.1:  The  geometry  of  the  2D  VXF  signal  propagation  problem.  The  properties 
of  the  ionosphere  and  the  Earth  are  assumed  to  be  constant  in  the  y  direction. 


The  customary-  means  for  entering  input  parameters  into  the  LWPC  code  is  in 
terms  of  a  model  file  which  contains  all  of  the  required  parameters  for  different  com¬ 
ponents  of  the  LWPC.  This  model  file  is  generally  named  TTTRR.lwpc  where 
TTT  stands  for  the  three  letters  describing  the  transmitter,  and  RR,  stands  for  the 
two  letters  describing  the  receiver.  After  the  input  file  is  prepared,  a  script  named 
modeL2d  is  used  to  run  the  different  FORTRAN  programs  described  below  in  the 
required  order.  Figure  2.2  shows  the  flow  of  this  process. 

Some  of  the  parameters  listed  in  the  TTTRR.Iw'-pc  file  and  used  by  different 
components  of  LWPC  are  listed  below  : 


/.  frequency  in  kHz  (default  =  23.4), 

P.  radiated  power  in  kW  (default  =  1), 

Ticmq-Tiat,  transmitter  coordinates  in  degrees  west  and  north  (default  =158.15/21.41), 
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Ri<mg,Riat,  receiver  coordinates  in  degrees  west  and  north  (default  =  0/0), 

dp  distance  increment  in  km  (default  =  5), 

Pmax,  maximum  distance  in  km  (default  =  10000), 

year, month, day,  for  terminator  calculations  (default  =  1984/03/15), 

gmt  Greenwich  time  in  hours  (default  =  -24), 

nrseg  number  of  segments  in  terminator  boundary  (default  =  5), 

ranger,  rangd,  range  of  real  and  imaginary  parts  of  the  eigenangles  to  be  found 
(default  =  60.0/90.0,  0.0/-7.0), 

Omax  maximum  attenuation  rate  in  dB/Mm  used  by  MODESRCH/MODEFNDR  to 
identify  a  finite  number  of  most  significant  modes,  (default  =  200), 

X  terminator  boundary  in  deg  of  solar  zenith  angle  (default  =  90), 

Pcap  polar  cap  boundary  in  deg  of  magnetic  dip  angle  (default  =  70), 

Rcomp,  received  electric  field  component,  =1  if  2,  =2  if  y  (default  =  1), 

Rait,  Tait  receiving  and  transmitting  platform  altitude  in  km  (default  =  0.0), 

y,  inclination  of  the  transmitting  antenna  with  respect  to  the  vertical  in  degrees 
(default  =  0.0), 

9,  orientation  of  the  transmitting  antenna’s  projection  in  the  horizontal  plane  with 
respect  to  the  direction  of  propagation  in  degrees  (default  =  0.0),  and 

Topht,  height  of  the  bottom  of  the  ionosphere  in  km  (default  =  90). 
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Figure  2.2:  Flow  diagram  for  the  LWPC  code  for  2D  modeling  of  VLF  signal  prop¬ 
agation  in  the  earth-ionosphere  waveguide.  TTTRR.lwpc  is  the  input  model  file 
which  contains  the  required  parameters  for  the  three  FORTRAN  programs,  namely 
PRESEG,  MODEFINDER  and  FASTMC,  described  in  the  text. 

2.1.1  PRESEG 

PRESEG  is  a  FORTRAN  program  which  uses  the  TTTRR.lwpc  model  file  and 
segments  the  propagation  path  between  the  transmitter  and  the  receiver  based  on 
ionosphere,  ground,  and  magnetic  field  changes  [Pappert  and  Snyder,  1972a;  Wait  and 
Spies,  1968].  It  creates  output  files  mf_TTTRR.com  and  fmc_TTTRR.com  which 
contain  this  segmentation  information  and  which  are  used  by  MODEFINDER  and 
FASTMC  respectively.  PRESEG  segments  the  path  according  to  the  directions  given 
in  the  TTTRR.lwpc  file  as  well  as  automatic  segmentation  rules,  such  as  the  known 
(as  derived  from  a  previously  input  conductivity  map  of  the  Earth)  conductivity 
changes  of  the  Earth  [Hauser  et  al,  1969]  as  the  signal  propagates  along  the  great 
circle  path,  crossing  boundaries  such  as  those  between  ground  and  sea,  rocky  versus 
marshy  soil,  as  well  as  different  ionospheric  regions,  such  as  low  or  mid  latitudes, 
auroral  regions,  etc.. 
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2.1.2  MODEFINDER 

The  solutions  to  the  mode  equation  (2.8)  in  general  cannot  be  expressed  in  closed 
form  for  realistic  ionospheric  conditions.  For  this  purpose,  a  computer  program  called 
MODEFINDER  is  employed  [Morfitt  and  Shellman,  1976;  Ferguson  and  Snyder, 
1987].  This  program  was  developed  initially  by  the  Naval  Ocean  Systems  Center 
in  order  to  numerically  determine  the  values  of  the  sine  of  the  complex  eigenangle, 
5n,  as  defined  in  connection  with  equation  (1.9).  It  is  designed  to  obtain  these  mode 
constants  rapidly  and  accurately.  The  MODEFINDER  program  obtains  the  full-wave 
solutions  of  the  modal  equations  for  —  Fq  in  terms  of  the  eigenangles. 

Two  other  input  parameters  that  are  included  in  the  TTTRR.lwpc  model  file 
are  the  altitude  profiles  of  the  electron  density  and  collision  frequency.  The  MOD¬ 
EFINDER  program  uses  this  information  and  first  modifies  the  density  profile  in 
order  to  eliminate  those  heights  for  which  the  electron  densities  have  negligible  ef¬ 
fects  on  the  modal  solutions.  The  modal  solutions  of  equation  (1.12)  are  determined 
within  a  bounded  region  of  the  complex  6  plane  determined  by  ranger,  rangd,  and 
ftmax-  In  the  process  of  determining  the  reflection  coefficients  of  the  ionosphere,  MOD¬ 
EFINDER  calculates  an  effective  reflection  height.  From  the  resultant  values  of  the 
reflection  coefficient  matrix  elements  and  the  modal  eigenangle  solutions,  the  attenu¬ 
ation  rate,  phase  velocity,  initial  excitation,  height  gain  factors  and  the  polarization 
mixing  ratios  for  each  mode  are  determined  [Morfitt  and  Shellman,  1976]. 

2.1.3  FASTMC 

In  order  to  calculate  the  mode  conversion  coefficient  matrices  at  the  different  slabs 
created  by  PRESEG,  a  FORTRAN  program  called  FASTMC  is  used  [Pappert  and 
Shockey,  1972;  Pappert  and  Snyder,  1972b].  FASTMC  is  a  simplified  version  of  a 
more  fundamental  FORTRAN  program  called  FULLMC.  FASTMC  uses  a  number 
of  approximations  in  order  to  reduce  the  execution  time  of  FULLMC  [Ferguson  and 
Snyder,  1980,  Poulsen,  1991]  and  is  exclusively  used  in  this  thesis  for  the  calculation 
of  mode  conversion  matrices. 
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2.1.4  Example  2D  simulation 

Since  the  LWPC  code  is  capable  of  accounting  for  mode  conversion  effects  at  the 
interface  between  two  slabs  having  one  or  more  different  geophysical  parameters,  we 
can  model  ionospheric  disturbances  which  lie  on  the  great  circle  propagation  path 
of  the  VLF  signal  propagating  from  a  transmitter  to  a  receiver,  as  long  as  such 
disturbances  can  be  represented  as  infinite  strips  perpendicular  to  the  path  and  with 
no  variation  in  the  y  direction.  Thus,  it  is  possible  to  study  the  effects  of  a  disturbance 
on  the  VLF  amplitude  and  phase  as  a  function,  for  example,  of  the  width  of  the 
disturbance  along  the  x  axis.  Conversely,  it  should  in  principle  be  possible  to  solve 
the  inverse  problem,  i.e.,  to  determine  the  width  (along  the  path)  of  a  disturbance 
from  measurements  of  the  associated  amplitude  and  phase.  In  Chapter  3,  we  employ 
such  an  inversion  method  to  determine  the  ambient  electron  density  profile. 

As  an  example  of  the  use  of  LWPC  for  2D  modeling  of  VLF  propagation,  we  con¬ 
sider  the  propagation  path  between  the  NLK  (24.8  kHz)  transmitter  in  Washington 
(121. 9°W,  48.20°N)  to  Fort  Yukon  (FY),  Alaska  (145.2'’W,  66.6°N).  We  also  assume 
the  hypothetical  presence  of  an  ionospheric  disturbance  along  the  propagation  path. 
The  NLKFY.lwpc  model  file  used  for  this  case  is: 

/export /home/demirkol/Pr ogr ams /nlkf y 
2 

/export /home/demirkol/profiles/ densprof s/electr .prof ile2 
/export /hpme/demirkol/prof  iles/densprofs/sp_relative_profs/d3_ll  .dat 
0 

NLKF 

NLKFY 

fedatum 

freq=24.8  trlat=48.200  trlong=121 . 910 
rclat=66.560  rclong=145.218 
max_alt=110.0 
iboundO  =  1400  1900 
inspht3rpe0=(l,0)  (2,0)  (1,0) 

&end 
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The  first  line  indicates  the  path  where  all  of  the  simulation  files  should  be  read  from 
and  written  to.  The  second  line  is  the  number  of  different  electron  density  profiles 
(in  this  case  2)  that  will  be  used  in  the  propagation  path  of  the  VLF  signal.  Next 
the  locations  of  these  tabular  electron  density  profile  files  are  given.  Alternatively,  it 
is  possible  to  use  an  exponential  profile  similar  to  that  in  equation  (1.2).  The 
number  0  on  the  next  line  corresponds  to  the  number  of  collision  firequency  profile  to 
be  used  in  the  calculations,  in  this  case  meaning  that  an  ambient  profile  should  be 
employed  in  this  specific  simulation.  The  next  few  lines  specify  the  frequency  of  the 
VLF  signal  as  well  as  transmitter  and  receiver  locations.  IBOUNDO  corresponds  to 
the  great  circle  distances  (from  the  transmitter)  in  km  where  a  user  forced  boundary 
exists,  in  this  case  representing  the  beginning  and  end  of  a  disturbance  between 
1400  to  1900  km.  INSPHTYPEO  corresponds  to  different  combinations  of  (electron 
density,  collision  frequency)  profiles  that  will  be  used  in  each  of  these  segments.  In  this 
case,  the  profile  combination  (1,0)  is  used  in  the  region  between  the  transmitter  and 
the  1400*’’  km  along  the  propagation  path,  consisting  of  the  first  (1)  electron  density 
profile  listed  in  the  model  file  and  the  ambient  (0)  collision  frequency  profile.  Between 
the  1400*’*  and  lOOO***  km  the  second  (2)  electron  density  profile  is  used  together  with 
the  ambient  (0)  collision  frequency  profile  corresponding  to  (2,0).  After  the  second 
boundary,  the  combination  is  again  (1,0).  The  geometry  of  the  problem  is  shown  in 
Figure  2.3. 

The  eigenangle  solutions  in  the  first  segment  of  the  propagation  path  are  shown  in 
Figure  2.4.  The  field  amplitude  and  phase  of  the  total  VLF  signal  at  different  points 
along  the  GCP  between  NLK  and  FY  are  shown  in  Figure  2.5.  The  characteristic 
deep  nulls  along  the  propagation  path  occurs  as  a  result  of  destructive  interference 
between  different  waveguide  modes,  each  propagating  with  slightly  different  phase 
and  velocity  and  different  attenuation  rates.  At  large  distances,  these  large  variations 
in  the  total  field  intensity  are  reduced,  since  most  of  the  modes  are  attenuated  out, 
and  the  signal  is  constituted  by  only  a  few  modes. 
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Ionosphere 


2.2  3D  Numerical  Modeling 

A  2D  model  is  necessarily  limited  in  application  since  most  realistic  ionospheric  dis¬ 
turbances  encountered  in  practice  have  a  three  dimensional  structure.  Examples  of 
3D  ionospheric  disturbances  are  those  produced  by  lightning  discharges  [Rodriguez  et 
ai,  1992]  and  powerful  HF  ionospheric  heaters.  A  3D  model  is  needed  for  accurate 
description  of  VLF  propagation  and  scattering  in  the  presence  of  such  disturbances. 
Such  a  3D  model  was  developed  at  Stanford  University  [Poulsen,  1991]  and  have  been 
applied  for  interpretation  of  VLF  signal  amplitude  and  phase  changes  observed  in  as¬ 
sociation  with  lightning  discharges  [Poulsen  et  al.  1993;  Lev-Tov  et  al,  1995;  Inan  et 
al,  1996].  The  geometry  of  the  3D  problem  is  shown  in  Figure  2.6.  The  3D  nature 
of  the  model  allows  for  ionospheric  parameter  variations  both  in  the  direction  of  the 
GCP  {x  direction)  and  the  direction  transverse  to  it  {y  direction) . 

In  the  3D  case  the  total  electric  field  a  distance  d  from  the  transmitter  as  given 
for  2D  by  equation  (1.9)  changes  to  equation  (2.1)  given  below  [Wait,  1964]  ; 
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Figure  2.4:  Eigenangle  solutions  found  by  MODEFINDER  for  the  NLK  signal  prop¬ 
agating  to  FY  for  the  first  segment  of  the  path.  Each  point  corresponds  to  a  different 
mode  as  indicated. 
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Figure  2.5:  The  amplitude  and  phase  of  the  VLF  signal  changes  as  it  propagates  in 
the  earth-ionosphere  waveguide.  The  solid  lines  correspond  to  the  forced  and  natural 
boundaries  on  the  propagation  path  as  segmented  by  PRESEG. 
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Figure  2.6:  Illustration  of  the  3D  modeling  of  the  VLF  wave  propagation  in  the 
presence  of  a  localized  disturbance.  The  main  leg  corresponds  to  the  GCP  between 
the  transmitter  and  the  receiver.  The  signal  at  the  receiver  consists  of  a  vector  sum 
of  the  direct  signal  along  the  Main  Path  and  the  scattered  signal  arrving  from  Leg  2. 


■^total  (d)  — 


d 

Re 


isin  (d/i?E)| 


exp 


-iko  j 


Sn  {x',  y')  ds 


(2.1) 


where  the  integration  contour  is  along  the  path  of  minimum  total  phase  between  the 
transmitter  and  the  observation  point  located  at  distance  d.  However  this  equation 
is  only  an  approximation,  assuming  slow  variation  of  the  properties  in  the  x  and  y 
direction.  A  more  detaiiled  analysis  of  this  problem  is  given  elsewhere  [Poulsen,  1991]. 

Based  on  the  formulation  described  in  Poulsen  [1991]  and  Poulsen  et  al.  [1993], 
the  3D  VLF  code  consists  of  the  LWPC  code  combined  with  two  other  programs  called 
SCATTER  and  SUMOUT,  in  order  to  simulate  the  propagation  of  VLF  signals  in 
the  presence  of  localized  lower  ionospheric  disturbances.  A  schematic  flow  graph  of 
the  3D  simulation  code  is  given  in  Figure  2.7,  where  SCATTER  and  SUMOUT  are 
briefly  described  below. 
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Figure  2.7:  Flow  of  the  LWPC  code  for  3D  simulation  of  VLF  signal  propagation  in 
the  earth-ionosphere  waveguide.  In  this  model,  VLF  signal  propagation  is  simulated 
on  two  different  legs  denoted  as  Main  and  Leg  1  in  Figure  2.6.  The  scattered  field  is 
calculated  by  the  FORTRAN  program  SCATTER  and  another  simulation  is  carried 
out  for  the  propagation  of  this  scattered  field  to  the  receiver,  over  the  path  denoted 
as  Leg  2.  Finally  SUMOUT  calculates  the  total  electric  field  that  is  received. 
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2.2.1  SCATTER 

The  program  SCATTER  calculates  the  theoretical  value  of  the  change  in  the  field 
measured  at  the  receiver  due  to  a  cylindrically  symmetric  perturbation  normalized 
bj  the  theoretical  value  of  the  field  that  would  have  been  received  had  there  been  no 
perturbation.  SCATTER  gives  both  the  real  and  imaginary  parts  of  this  change  in 
the  field.  A  few  of  the  important  parameters  in  this  program  are  : 

6o  the  background  complex  eigenangle  value, 

dpeak  the  eigenangle  values  at  the  center,  or  peak  of  the  disturbance, 

spotsize,  the  size  of  the  disturbed  region  in  the  ionosphere, 

XQ,yo  the  location  of  the  spot  with  respect  to  the  great  circle  path  of  VLF  signal 
propagation. 


This  program  uses  a  specified  lateral  profile  (e.g.,  a  Gaussian  distribution)  as  the 
horizontal  roll-off  function  for  the  change  in  S  value  of  the  perturbation  back  to  the 
background  value  as  a  function  of  distance  from  the  center  of  the  perturbation. 

2.2.2  SUMOUT 

The  total  field  measured  at  the  receiver  is  a  sum  of  the  field  from  the  main  leg,  which 
is  the  direct  propagation  of  the  signal  along  the  great  circle  path  and  the  scattered 
field  from  the  disturbance  arriving  over  Leg  2.  SUMOUT  implements  the  vectorial 
addition  of  two  fields. 

2.3  Analysis  of  a  Single  Mode  Case 

The  subject  of  this  dissertation  is  to  utilize  the  tools  just  described  in  order  to  measure 
the  electron  collision  frequency  profile  within  a  heated  ionospheric  patch  by  exploit- 
ing  the  sensitive  dependence  of  the  VLF  signal  amplitude  and  phase  on  ionospheric 
properties  such  as  electron  density  and  temperature.  In  this  section,  we  illustrate  the 
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determination  of  the  collision  frequency  profile  for  the  simplest  possible  case  in  which 
the  VLF  signal  is  constituted  by  a  single  waveguide  mode.  This  case  can  be  treated 
analytically  and  its  analysis  illustrates  the  difficulty  involved  in  the  more  general 
multimode  case. 

The  individual  modes  may  be  identified  in  terms  of  the  number  of  minima  and 
maxima  that  the  field  quantities  exhibit  as  a  function  of  altitude  (Figure  1.4-b)  or 
in  terms  of  the  discrete  eigenangles  for  which  the  mode  fields  can  be  coherently 
decomposed  into  a  sum  of  two  uniform  plane  waves  which  are  in  phase.  The  complex 
eigenangles  may  be  obtained  through  a  full  wave  solution  of  Maxwell’s  equations 
which  5delds  a  set  of  reflection  coefficients  satisfying  the  condition: 


r{9)  exp  {-2ikhsmd)  =  1 


(2.2) 


which  is  simply  obtained  from  equation  (1.12)  by  settinj^i 


Fo 


In  general  R^^  is  a  matrix  because  of  the  anisotropic  characteristic  of  the  iono¬ 
spheric  plasma.  However  during  daytime  when  the  VLF  reflection  height  approaches 
70  km  in  altitude,  the  effects  of  the  Earth’s  magnetic  field  are  reduced  and  the  iono¬ 
sphere  at  this  altitude  can  be  considered  isotropic.  In  this  case  is  a  scalar. 

In  general  the  reflection  coefficient  R{9)  is  the  ratio,  at  the  boundary,  of  some  field 
component  Fr  in  the  reflected  wave  to  the  same  field  component  Fi  in  the  incident 
wave.  Fr  and  Ft  are  both  measured  very  close  to  the  boundary  either  at  the  atmo¬ 
sphere  or  at  the  ground.  For  the  case  of  sharp  boundaries,  these  reflection  coefficients 
are  called  Fresnel  reflection  coefficients.  For  electromagnetic  waves  linearly  polarized 
with  the  electric  field  vector  in  the  plane  of  incidence  (the  x  —  z  plane)  the  reflection 
coefficient  at  the  ionospheric  boundary  is  given  as: 


R  = 


n  sin  ^  —  sin  ^ 
n  sin  0  -f-  sin  ^ 


(2.3) 


where 


n,  is  the  refractive  index  as  calculated  by  equation  (4.4), 

6  is  the  angle  between  the  wave  normal  in  the  x  —  z  plane  and  the  boundary. 
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^  is  the  angle  between  the  transmitted  wave  normal  and  the  boundary-. 
6  and  ^  are  related  to  each  other  by  Snell’s  law  : 


cos6  =  ncos'^  (2.4) 

For  electromagnetic  waves  linearly  polarized  with  the  electric  field  vector  perpen¬ 
dicular  to  the  plane  of  incidence  (i.e.,  the  x  —  z  plane)  the  reflection  coefficient  at  the 
ionospheric  boundary  is  given  as; 


_  sin  0  —  n  sin 
sin^-fnsin^' 

A  better  approximation  should  be  used  for  the  reflection  coefficient  from  the  iono¬ 
sphere  where  the  boundary  is  not  sharp.  In  this  case  a  stratified  medium  is  assumed 
where  the  refractive  index  n{z)  is  a  function  of  altitude,  equivalent  to  a  number  of  dis¬ 
crete  thin  layers  in  each  of  which  the  refractive  index  is  constant.  Snell’s  law  applies 
at  each  of  the  boundaries  : 


rife  cos  '^k  =  nfc+i  cos  ^*,+1  (2.6) 

In  this  case  it  can  be  shown  that  the  correct  formula  for  the  reflection  coefficient 
is  given  as  [Budden,  1961;  p.  117]: 


where 

q  =  \/rP  —  cos 


(2.7) 


2o  is  the  height  where  g  =  0  and  thus  at  which  the  wave  is  reflected. 

Under  night  time  conditions,  the  reflection  height  rises  to  approximately  90  km. 
where  the  Earth’s  magnetic  field  creates  a  strong  anisotropy  which  leads  to  polariza¬ 
tion  coupling.  In  this  case  the  reflection  coefficients  and  R‘  must  be  regarded  as 
matrices  and  written  in  the  form  R®  and  R'  where 
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±Rmd 
,,R±iO)d  ±R±{0)d 

trix  looking  up  in  the  ionosphere  from  height  d  and 


is  the  complex  ionospheric  reflection  coefficient  ma- 


.RK0)d  0 

0  ±Rl{e)d 

height  d  towards  the  ground 


is  the  complex  reflection  matrix  looking  down  from 


The  notation  ||  for  R's  denotes  polarization  parallel  to  the  plane  of  incidence  while 
the  notation  ±  denotes  polarization  perpendicular  to  the  plane  of  incidence.  The  first 
subscript  on  the  R's  refers  to  the  polarization  of  the  incident  wave  while  the  second 
applies  to  the  polarization  of  the  reflected  wave.  These  reflection  coefficients  are 
calculated  and  used  in  LWPC  simulations. 

In  this  case  the  modal  equation  given  in  (2.2)  converts  to  a  determinantal  equation 
given  by: 


=  0 


(2.8) 


The  total  electric  field  at  a  great  circle  path  (GCP)  distance  d  from  a  transmitter 
has  the  form 


E/total{,d^  — 


1 

>/|sin(d/RE)| 


OO 

exp  [-ffco  sin  {6n)  d] 

n 


(2.9) 


and  is  the  electric  field  that  would  be  measured  along  a  direct  path  between  a  trans¬ 
mitter  and  receiver.  In  the  presence  of  an  HF  heater  induced  disturbance  of  the  D 
region  of  the  ionosphere,  there  is  the  additional  field  scattered  from  the  disturbed  re¬ 
gion  which  causes  the  amplitude  and  phase  changes  in  the  total  received  VLF  signal 
to  be  different  than  that  given  by  equation  (2.9). 

We  now  consider  the  hypothetical  case  in  which  only  a  single  waveguide  mode 
(say  mode  n=l)  is  excited,  so  that  n=l  and  A-^’  =1.  Assuming  that  the  ambient 
electron  density  profile  and  collision  frequency  profile  are  known  the  field  amplitude 
and  the  eigenangle  corresponding  to  each  mode  can  be  calculated  using  equation  (2.2) 
and  (4.4).  Thus  the  direct  field  measured  at  a  receiver  station  in  the  absence  of  the 
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disturbance  is  given  by: 


where 


^°{^R,yR) 


exp  [— ifco  sin  (0°)  d] 

v/|sin(d/fiE)l 


(2.10) 


ko  =  w/c, 

6°  is  the  eigenangle  corresponding  to  the  dominant  single  mode  propagating  along 
the  direct  path  between  the  transmitter  and  the  receiver,  sin(^°)  is  generally 
simply  written  as  S°. 

d  =  y/ {xr  -  XtY  +  (j/fi  -  vtY, 

{xR,yR)  and  {xt^vt)  are  the  receiver  and  transmitter  locations  respectively, 
is  the  Earth  radius. 

Using  a  single  mode  scattering  formulation  [Poulsen  et  al,  1990]  the  scattered 
field  at  the  receiver  is  given  by; 


vr)  =  f{xR,  yR)  +  f{xR,  yR)  (2.ii) 

'  ^ . ^  'V'  '  ' 

direct  field  scattered  field 

=  e°{xR,  yR)  -  ^  [ j lS‘^{x' ,y')  -  (5'°)^]  ■  e(^^- 

^  field  inside  disturbance 

H^o\koS°R!)dx'dy' 

where 

e  is  the  modal  wave  field  (i.e.,  the  total  field  seen  at  {xr,  yR)  in  the  presence  of  some 
disturbance  in  the  waveguide) , 

e°  is  the  unperturbed  modal  wave  field,  also  called  the  ’direct  field’, 
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e®  is  the  secondary  or  scattered  field  (i.e.,  the  field  seen  at  {xfi,yR)  due  to  the 
disturbance), 

S°  is  the  .  ambient  value  of  S  in  the  absence  of  any  disturbance,  (in  the  single  mode 
case  this  is  a  constant  independent  of  x  and  y), 

AS^  the  difference  between  the  square  of  the  sine  of  the  perturbed  eigenangle  with 
respect  to  the  square  of  the  sine  of  the  ambient  eigenangle, 

P  is  the  region  of  integration  which  extends  over  that  portion  of  the  x  -  y  plane 
that  encompasses  the  disturbance, 

is  the  Hankel  function  of  second  kind  of  order  zero,  and 

R'  =  yJ{xR  -  x'Y  +  {ur  -  y'Y  is  the  distance  from  each  integration  point  (x/,  y/) 
within  the  disturbance  to  the  observation  point  (xH,yj?). 


The  amplitude  change  ratio  measured  at  the  receiver  station  is  then  simply  given 


by: 


-Ar  = 


\e°  (xR,yR)  +  e®  {xR,yR)\ 


\e°  {xR,yR)\ 

and  the  phase  change  measured  at  the  receiver  station  is  given  by: 


(2.12) 


A(j)  =  [e°  {xr,  yR)  +  e®  (xr,  yp)]  -  Ze°(xR,  yn)  (2.13) 

Since  the  value  e°{xR,yR)  can  be  calculated  with  (2.10)  for  a  specific 
transmitter-receiver  configuration  when  the  ambient  electron  density  and  collision 
frequency  profile  is  known,  it  is  in  principle  possible  to  find  the  complex  valued 
e*(xR,  yR)  by  using  amplitude  and  phase  changes  in  the  VLF  signal.  We  now  demon¬ 
strate  how  this  could  be  done  using  3  different  VLF  measurements,  as  described  in 
Chapter  1.  For  this  purpose  we  continue  to  make  the  hypothetical  assumption  that 
all  VLF  signals  (direct  and  scattered)  consist  only  of  one  waveguide  mode,  and  that 
we  have  in  hand  three  separate  measurements  of  (xr,  yR)  at  three  different  frequen¬ 
cies,  which  may  well  have  to  be  acquired  at  three  diflferent  locations  (see  Figure  4.3), 
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as  is  necessary  in  practice.  Since  the  three  calculated  scattered  fields  e^XR^yn)  for 
each  of  the  receiver  sites  are  complex,  this  information  can  be  used  to  determine  the 
collision  frequency  profile  with  6  parameters  or  with  6  degrees  of  freedom. 

There  are  mainly  two  steps  in  this  inversion  of  scattered  fields  to  determine  the 
collision  frequency  profile.  The  first  step  is  to  find  the  eigenangle  of  the  mode  corre¬ 
sponding  to  the  propagating  scattered  signal.  For  this  purpose,  we  make  use  of  the 
fact  [Poulsen  et  ai,  1993]  that: 


e^i^R^VR)  =  f[S{x\y')  ,P]  (2.14) 

=  -X  //  -  iSn"]e{x’,y')  H^^\koS°R')dx’dy'  (2.15) 

P 

where  e*  {xr,  yR)  is  a  function  of  the  unknown  S  (x',  y')  which  is  the  sine  of  the  eige¬ 
nangle  of  the  scattered  mode  and  P  corresponds  to  the  region  disturbed  by  the  heater. 
Equation  (??)  can  be  simplified  by  assuming  a  cylindrically  symmetric  disturbance 
with  a  center  at  (x^,  yj)  and  having  a  radius  of  ro  and 


S{x', y')  =  exp  ^  (sin -  sin 0°)  -I-  sin (2. 16) 


where  we  assume  a  Gaussian  form  (in  the  lateral  direction)  for  the  disturbance  which 
generates  a  scattered  signal  mode  with  an  eigenangle  value  of  6^. 

In  this  case,  e®(x/?,  yj? )  is  of  the  form  A  sin  9^ -B  sin  9°  where  A  and  B  are  constants 
that  can  be  found  by  doing  the  integration  given  in  (2.14).  Thus  sin^®  = 
sin  ^  {1/A  [e*  {xR.yR)  -I-  Bsin  (0°)]}  and  it  is  possible  to  directly  find  the  eigenangle 
9^  of  the  scattered  mode  for  three  different  VLF  signals. 

The  second  step  is  to  find  the  vector  u  =  (^'i,r'2)r'3,^'4,t'5,t'6)  >  the  vector  con¬ 
taining  the  six  values  of  the  disturbed  collision  frequency  profile  u{z)  at  six  different 
altitudes  Xj,  satisfying  the  model  equation  (2.8)  with  for  three  different  VLF 

signal  frequencies.  However  it  is  simpler  to  express  the  modal  equation  in  the  form  : 
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h{i;,e^)  =  o  (2.17) 

where  h  is  a.  function  of  different  parameters  of  which  only  P  is  unknown.  In  this 
form  it  is  easier  to  visualize  the  Newton-Raphson  algorithm  used  for  the  inversion 
and  described  below:. 

1)  Start  with  an  initial  estimation  of  > ^3 > ^4!  ^5 > ^e)  • 

2)  For  the  m*^  iteration,  let 


Sh  h{i7^  +  Sup,  01)  -h{V^-  Si^,  d]) 

Si/i  2A^'i 

where 

6ir  =  ociU..,  ^  ,..,0 

\  position 

j  =  1,2,3  corresponding  to  three  different  frequencies, 
i  =  1, 2, ..,  6  for  six  different  collision  frequency  values. 


(2.18) 


aj  is  a  constant  which  will  be  referred  as  step  size. 

3)  From  among  the  twelve  different  values  of  h  {p^  ±  6f)  for  i  =  1, 2,  ..,6  choose 

3 

i  such  that  s  =  ^i)  1  is  minimized.  Let  this  index  be  k. 

j=i 

4)  Update  the  value  of  +  Au,  i.e.,  either  increment  or  decrement  the  real 

or  the  imaginary  part  of  by  or  AOi  by: 


/ 


\ 


Ai7  = 


0,..,  \-.l-2/ai,..,0 


(2.19) 
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5)  Go  to  step  2  and  continue  until  e  <  Sdesired- 

It  is  apparent  from  this  single  mode  analysis  that  it  is  in  principle  possible  to  find 
the  collision  frequency  profile  within  a  disturbed  region  approximately  by  processing 
a  set  of  observed  VLF  amplitude  and  phase  changes.  In  practice,  VLF  signals  propa¬ 
gating  in  the  earth-ionosphere  waveguide  are  constituted  of  a  summation  of  multiple 
waveguide  modes,  so  that  the  total  field  at  the  receiver  in  the  absence  of  a  disturbance 
is  given  by  equation  (2.9),  and  the  eigenangles  ^  and  excitation  constants  A^^  must 
be  calculated  numerically  using  LWPC  as  described  in  Section  2.1.  The  scattered  sig¬ 
nal  at  the  receiver  must  also  be  calculated  numerically  using  the  3D  code  described 
in  Section  2.2.  The  general  algorithm  used  in  such  a  case  for  the  inversion  of  VLF 
measurements  to  determine  the  collision  frequency  profile  is  described  in  Section  4.2. 


Chapter  3 

Ambient  Electron  Density  Profile 

The  amplitude  and  phase  of  subionospheric  VLF  signals  are  sensitively  dependent  on 
the  altitude  profile  of  electron  density  along  the  propagation  path.  In  this  chapter 
we  pursue  the  possibility  of  utilizing  this  dependence  for  the  determination  of  the 
ambient  (i.e.,  prior  to  ionospheric  modification  via  HF  heating)  electron  density  from 
measured  values  of  VLF  amplitude  and  phase.  Knowledge  of  the  ambient  electron 
density  profile  greatly  simplifies  the  determination  of  the  collision  frequency  profile 
within  a  disturbed  region,  as  described  in  Section  4.1. 

The  latitude  region  of  interest  in  this  dissertation  is  the  subauroral  region  in  which 
the  HAARP  facility  is  located,  and  the  ionospheric  conditions  at  these  latitudes  are 
strongly  influenced  by  relativistic  electron  precipitation.  It  has  recently  been  shown 
[Demirkol  et  al,  1999]  that  the  amplitude  of  subionospheric  VLF  signals  propagat¬ 
ing  in  these  regions  are  indeed  closely  associated  with  relativistic  electron  fluxes  as 
measured  on  satellites.  This  close  association,  and  the  fact  that  relativistic  electron 
enhancements  on  the  average  occur  in  a  well  defined  range  of  subauroral  latitudes, 
enables  us  to  interrelate  the  ambient  electron  density  profiles  at  different  points  along 
the  subauroral  traverse  of  a  VLF  path,  which  in  turn  facilitates  the  determination  of 
the  ambient  density. 

In  the  following  we  review  the  results  reported  in  Demirkol  et  a/.  [1999]  and  describe 
the  method  used  to  determine  the  ambient  electron  density  profile  from  VLF  data. 
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3. 1  Ionospheric  Effects  of  Relativistic  Electron  En- 
hancement  Events 

The  Earth  s  outer  magnetosphere  is  often  populated  to  a  surprising  degree  by  rela¬ 
tivistic  electrons  [Paulikas  and  Blake,  1979].  Enhancements  in  the  relativistic  electron 
fluxes  may  be  an  important  source  of  energy  input  to  the  atmosphere.  Those  precip¬ 
itating  electrons  with  energies  >  1  MeV  can  penetrate  to  altitudes  as  low  as  50  km, 
affecting  the  atmospheric  chemistry  throughout  the  mesosphere  [Games  et  al,  1995]. 
Relativistic  electron  precipitation  events  are  also  believed  to  be  a  significant  source 
of  odd  nitrogen  in  the  middle  atmosphere,  possibly  affecting  ozone  concentrations  in 
some  regions  of  the  atmosphere  [Callis  et  ai,  1991].  Relativistic  electron  precipitation 
events  are  associated  with  magnetospheric  activity  and  may  appear  more  firequentiy 
near  a  solar  minimum  then  solar  maximum  [Baker  et  al,  1986].  These  events  are 
strongest  at  subauroral  (4.5  <L<7)  latitudes  (See  Figure  3.1). 

The  enhancements  of  energetic  particle  fluxes  within  and  near  the  local  loss  cone 
are  documented  in  data  from  low  altitude  satellites  such  as  SAMPEX  [Baker  et  al. 
1986]  and  UARS  [Gaines  et  al,  1995],  while  the  relativistic  electron  population  at 
geosynchronous  orbit  is  measured  on  GOES-7  and  GOES-8.  The  particle  flux  as  mea¬ 
sured  on  these  satellites  exhibit  the  well  known  relatively  regular  27-day  periodicity 
with  typical  rises  on  a  2  to  3-day  time  scale  and  decays  on  a  3  to  4-day  scale  [Baker 
et  al,  1986].  (See  Figure  3.2). 

As  mentioned  previously,  VLF  sounding  of  the  lower  ionosphere  (i.e.,  the  mea¬ 
surement  of  the  amplitude  and  phase  of  subionospherically  propagating  VLF  signals) 
is  a  sensitive  tool  for  the  detection  of  ionospheric  conductivity  changes  due  to  changes 
m  electron  density  and/or  temperature,  especially  at  altitudes  below  90  km  [Sechrist, 
1974].  Some  of  the  early  work  on  relativistic  electron  precipitation  events  has  in¬ 
deed  relied  on  subionospheric  VLF  measurements  [Thorne  and  Larsen,  1976].  Thus 
VLF  measurements  can  be  used  to  quantitatively  assess  the  degree  to  which  rel¬ 
ativistic  electron  enhancements  observed  at  satellite  altitudes  are  accompanied  by 
enhanced  precipitation  into  the  ionosphere.  For  this  purpose,  we  search  for  the  iono¬ 
spheric  signatures  of  relativistic  electron  precipitation  by  interpreting  the  observed 
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Figure  3.1:  Configuration  of  NLK-FY  path  with  resp«:t  to  SAMPEX  measurements 
of  the  relatmstic  precipitation  zone.  The  NLK-FY  path  is  chosen  since  the  strongest 
fiux  levels  are  at  the  subauroral  latitudes  traversed  by  this  path. 
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Figure  3.2:  The  particle  flux  measured  on  the  GOES-7  and  SAMPEX  satellites  ex¬ 
hibits  the  well  known  relatively  regular  27-day  periodicity  with  typical  rises  on  a  2  to 
3-day  time  scale  and  decays  on  a  3  to  4-day  scale. 
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VLF  amplitude  variations  in  the  light  of  theoretical  models  of  VLF  subionospheric 
wave  propagation  described  in  Chapter  2.  Our  results  indicate  that  the  nighttime  D 
region  is  indeed  strongly  affected  by  this  precipitation,  with  the  nighttime  electron 
density  at  40-70  km  altitudes  clearly  exhibiting  the  27-day  cycle  associated  with  solar 
rotation. 

3.1.1  Description  of  data 

For  this  study,  we  focus  our  attention  on  the  two  month  period  of  October  and 
November  1992  (for  which  VLF  data  from  Fort  Yukon  is  continuously  available)  and 
utilize  three  different  data  sets,  namely  (i)  VLF  data,  (ii)  SAMPEX  data,  and  (iii) 
GOES  7  data. 

The  VLF  data  consists  of  the  recorded  amplitude  and  phase  of  the  subionospheric 
signal  from  the  NLK  transmitter  (24.8  kHz)  in  Jim  Creek,  Washington  (121.91°  W, 
48.20°  N)  as  received  at  Fort  Yukon  (FY),  Alaska  (145.21°  W,  66.56°  N)  during  the 
period  Oct-Nov  1992.  Figure  3.3-a  shows  the  NLK-FY  great  circle  propagation  path  as 
well  as  lines  of  constant  geomagnetic  latitude.  Since  the  fluxes  of  relativistic  electrons 
exhibit  fluctuations  primarily  at  subauroral  latitudes  (4.5  <  L  <  7),  the  NLK-FY 
path  is  well  situated  for  monitoring  the  ionospheric  affects  of  these  relativistic  electron 
enhancements. 

Increases  in  the  D  region  electron  density  caused  by  the  high-energy  particle 
precipitation  increases  the  local  electrical  conductivity  and  perturbs  the  VLF  signal 
propagating  under  the  disturbed  ionosphere.  Figure  3.3-b  shows  a  schematic  descrip¬ 
tion  of  this  process. 

The  VLF  data  at  FY  during  the  Oct-Nov  1992  period  was  typically  recorded 
during  the  period  0000  to  1200  UT.  The  signal  amplitude  in  a  300  Hz  band  centered 
at  the  transmitter  frequency  (24.8  kHz)  is  regularly  sampled  and  digitally  recorded  at 
a  resolution  of  100  Hz  (i.e.,  samples  taken  at  10  ms  intervals).  Since  the  VLF  signal 
amplitude  can  exhibit  significant  variation  over  short  time  scales,  for  example  in 
response  to  burst  precipitation  effects  [Cotton  and  Smith,  1991]  or  auroral  electrojet 
enhancements  [Cummer  et  ai,  1996],  studies  of  long  term  behavior  are  facilitated 
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Figure  3.3:  (a)  A  geographic  view  of  the  VLF  propagation  path  from  NLK  to  Fort 
Yukon,  Alaska.  The  NLK-FY  path  is  sitiiated  such  that  the  relativistic  precipi¬ 
tation  region  (shown  shaded)  covers  a  significant  portion  of  the  great  circle  prop¬ 
agation  path,  (b)  The  mechanism  of  the  particle  precipitation-VLF  interaction  is 
schematic^y  shown  above.  A  ground  VLF  transmitter  (T)  launches  a  signal  into 
the  earth-ionosphere  waveguide.  In  the  region  of  relativistic  particle  precipitation  the 
local  electron  density  is  increased  by  A  ATg  .causing  the  local  electrical  conductivity  to 
change,  “nie  wavegmde  signal  propagating  under  the  region  of  relativistic  electron 
precipitation  is  modified  in  response  to  the  conductivity  change,  allowing  the  ob¬ 
servation  of  this  conductivity  change  as  phase  and  amplitude  variations  in  the  VLF 
signal. 
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through  the  use  of  data  averaged  over  a  number  of  hours  on  each  day.  For  this 
purpose  we  chose  to  average  the  VLF  data  over  the  time  interval  0600  to  0900  UT 
each  day,  during  which  the  entire  NLK-FY  path  was  under  darkness  throughout  the 
study  period.  For  each  day,  this  simple  3  hour  average  of  the  signal  amplitude  is 
designated  as  the  ‘average’  signal  amplitude. 

The  Proton/Electron  Telescope  (PET)  [Cook  et  al,  1993]  on  SAMPEX  is  com¬ 
posed  of  an  array  of  silicon  soUd  state  detectors  that  identify  and  measure  the  kinetic 
energy  of  electrons  from  ~1  to  ~30  MeV  and  of  H  and  He  isotopes  from  ~20  to  ~80 
MeV/nuc.  The  SAMPEX  data  used  in  this  work  are  integral  fluxes  of  >  4  MeV 
electrons  measured  in  specific  passes  nearby  the  subionospheric  paths. 

The  GOES-7  Energetic  Particle  Sensor  (EPS)  measures  electrons  from  0.6  to 
greater  than  4.0  MeV,  protons  from  .8  to  500  MeV,  and  aJpha-particles  from  4  to 
500  MeV  [Goes  Handbook,  1994]  .  The  electron  measurements  are  made  via  solid 
state  surface  barrier  detectors  within  a  dome  subassembly.  The  data  used  in  this 
work  are  5-minute  averages  of  the  integral  fluxes  of  >  2  MeV  electrons. 

3.1.2  Ground  and  Satellite  Data  Comparisons 

The  ‘average’  VLF  amplitude  for  the  October  and  November  1992  period  is  compared 
in  Figure  3.4  with  the  corresponding  satellite  particle  data  measured  by  GOES-7.  The 
GOES-7  electron  flux  maxima  are  clearly  associated  with  the  amplitude  minima  in 
the  VLF  data.  A  27-day  cycle  is  clearly  apparent  in  the  VLF  data,  similar  to  that 
observed  in  GOES  data  for  the  relativistic  electron  enhancement  events,  with  the 
same  3-4  days  of  rise  and  fall  times.  In  association  with  the  two  large  electron  flux 
enhancements  which  peak  on  days  276  and  305,  the  VLF  signal  amplitude  exhibits 
changes  >  9  dB  with  rise  and  fall  times  of  a  few  days.  The  VLF  amplitude  minima 
are  delayed  by  about  two  days  with  respect  to  the  peaks  in  GOES  data;  possible 
reasons  for  this  delay  are  discussed  later.  The  two  smaller  peaks  in  the  GOES  data, 
which  are  an  order  of  magnitude  below  the  two  main  peaks,  are  not  associated  with 
strong  VLF  amplitude  minima. 

To  our  knowledge  this  is  the  first  observation  of  a  subionospheric  VLF  amplitude 
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Days  of  1992 
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Figure  3.4:  (a)  The  three-hour-averaged  amplitude  of  the  NLK  signal  (24.8  kHz) 
at  Fort  Yukon  (FY)  on  each  day.  (b)  Electron  precipitation  flux  as  measured  on 
GOES-7.  A  27-day  variation  is  apparent  in  each  data  set. 
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variation  exhibiting  the  same  27-day  cycle  as  the  relativistic  electron  enhancements 
events,  thus  indicating  that  the  nighttime  lower  ionospheric  electron  density  is  de- 
tectably  affected  by  the  solar  rotation.  We  note  that  since  we  use  3  hr  averages  of 
the  VLF  data,  any  measurement  noise  in  the  300  Hz  VLF  channel  is  significantly 
reduced,  and  the  daily  variations  shown  in  Figure  3.4-a  are  true  indications  of  day  to 
day  changes  in  the  lower  ionosphere,  principally  due  to  auroral  effects  and  long  term 
precipitation  associated  with  the  auroral  electrojet  [Cummer  et  al.,  1996].  The  27 
day  variation  is  clearly  the  dominant  effect,  imposed  on  top  of  these  other  variations. 
At  the  same  time  the  variability  of  the  nighttime  ionosphere  due  to  the  other  auroral 
effects  probably  accounts  for  the  lack  of  a  VLF  amplitude  minima  associated  with  the 
smaller  peak  in  the  GOES-7  data  near  day  292  (Figure  3.4).  Using  the  three-hour- 
averaged  amplitude,  it  appears  that  only  relativistic  electron  enhanced  flux  levels 
above  ~3xl0^  el/cm^-sr-s  produce  ionospheric  effects  that  stand  out  in  the  presence 
of  other  ionospheric  variations.  Figure  3.5  shows  the  VLF  signal  amplitude  for  the 
same  three  hour  period  for  two  days  when  the  flux  level  is  at  the  maximum  and 
minimum  levels. 

In  view  of  the  remarkable  clarity  by  which  the  relativistic  electron  enhancement 
peaks  are  manifested  in  the  daily  averaged  VLF  data,  it  is  interesting  to  consider 
whether  the  NLK-FY  and  GOES-7  data  sets  exhibit  correlation  on  shorter  time  scales. 
Figure  3.6shows  comparison  of  data  at  5  min  resolution  for  a  selected  period  (0600- 
0900  UT  on  October  25,  1992)  in  the  vicinity  of  one  of  the  peaks.  We  see  that  the 
VLF  intensity  increases  in  response  to  a  decrease  in  the  electron  flux.  In  general  a 
correlation  coefficient  of  magnitude  ~0.7  was  calculated  between  the  two  time  series. 

3.1.3  Model  Calculations 

The  coincident  occurrence  of  subionospheric  VLF  signal  changes  and  relativistic  elec¬ 
tron  enhancement  peaks  suggest  that  significant  enhanced  precipitation  accompanies 
the  enhancement  events.  In  order  to  determine  whether  the  observed  VLF  amplitude 
signatures  are  consistent  with  the  ionospheric  changes  expected  to  be  produced  by 
such  relativistic  electron  precipitation,  we  theoretically  model  the  propagation  of  the 
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Figure  3.5:  VLF  signal  amplitudes  on  days  300  and  307  shown  in  the  interval  0600- 
0900  UT.  The  difference  in  the  amplitude  for  these  two  days  is  >  10  dB.  VLF  signal 
amplitude  from  Day  307  corresponds  to  a  day  where  an  electron  precipitation  en¬ 
hancement  takes  place,  while  the  signal  from  Day  300  is  a  day  where  the  relativistic 
electron  flux  level  is  moderate. 
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«  E  >  2  MeV 


Figure  3.6:  5-ininute  average  GOES-7  data  appears  to  be  strongly  correlated  with 
VLF  amplitude.  Correlation  coefficients  up  to  ~  —0.7  are  calculated. 
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VLF  signal  in  the  earth-ionosphere  waveguide  along  the  great  circle  path  from  NLK 
to  FY.  For  this  purpose,  we  use  the  Long  Wave  Propagation  Capability  (LWPC) 
code  described  in  Section  2.1.  The  disturbance  along  the  NLK-FY  propagation  path 
is  modeled  as  a  segment  of  the  earth-ionosphere  waveguide  with  a  perturbed  elec¬ 
tron  density  profile  caused  by  relativistic  electron  precipitation.  Figure  3.7-a  and 

3.7- b  show  perturbed  electron  density  profiles  associated  with  diflferent  flux  levels  as 
determined  by  the  method  of  Gaines  et  al.  [1995].  Energetic  electron  data  from 
SAMPEX,  specifically  integral  flux  of  >  4  MeV  electrons,  which  was  recorded  during 
passes  when  the  satellite  was  closest  (in  longitude)  to  the  propagation  path,  is  used  in 
order  to  have  the  best  estimate  of  electron  flux  levels  in  the  ionospheric  region  which 
lies  above  the  VLF  path.  A  time  period  in  the  vicinity  of  the  peak  in  relativistic  flux 
on  day  305  was  chosen  for  this  purpose. 

There  are  two  notable  characteristic  of  the  flux  enhancements  shown  in  Figure 

3.8- a.  First,  we  note  that  the  largest  flux  level  changes  occur  between  5  <  L  <  7. 
Second,  we  note  a  significant  spatial  expansion  of  the  disturbed  region  following  day 
304  which  covers  an  increasingly  larger  segment  of  the  NLK-FY  path.  This  spatial 
expansion  is  the  likely  reason  for  the  ~2  day  delay  of  the  VLF  amplitude  minima  with 
respect  to  the  relativistic  enhancement  peak  which  was  noted  earlier  in  connection 
with  Figure  3.4.  Note  that  the  VLF  amplitude  change  is  proportional  to  the  electron 
enhancement  as  well  as  the  length  of  the  path  segment  affected  [Inan  and  Carpenter, 
1987]. 

For  the  VLF  propagation  model  calculations,  the  propagation  path  NLK-FY  was 
segmented  using  electron  density  profiles  associated  with  the  actual  flux  levels  shown 
in  Figure  3.8-a.  Figure  3.8-b  shows  flux  levels  used  in  our  model  as  derived  from  SAM¬ 
PEX  data  for  I-~6.  LWPC  code  calculations  were  carried  out  using  these  models  of 
the  ionospheric  disturbance.  The  results  shown  in  Figure  3.9  predict  ~7  dB  maximum 
amplitude  decrease,  which  compares  generally  well  with  the  observed  9  dB  signal  am¬ 
plitude  decrease.  The  somewhat  lower  calculated  amplitude  change  may  be  due  to 
the  fact  that  >  4  MeV  electron  flux  used  to  determine  the  associated  electron  den¬ 
sity  profiles,  underestimated  the  relativistic  electron  enhancements  which  generally 
involve  electrons  with  energies  >  0.5  MeV. 
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Figure  3.7:  (a)  Ionization  profiles  corresponding  to  the  different  levels  (1,  2,  3, 4,  5)  of 
relativistic  electron  precipitation  fluxes  as  shown  in  (b)  A  typical  relativistic  electron 
precipitation  enhancement  shown  here  rising  and  falling  in  9  days.  In  many  cases 
event  durations  are  as  much  as  10-15  days.  The  flux  levels  and  energy  spectra  of  the 
precipitation  was  taken  to  be  as  given  by  Gaines  et  al.  [1995],  based  on  measurements 
on  the  UARS  satellite. 


Flux  variation  on  consecutive  days 


Figure  3.8:  (a)  Integral  flux  of  >  4  MeV  electrons  measured  on  each  day  on  SAMPEX 
during  passes  nearest  (in  longitude)  to  the  NLK-FY  path.  The  largest  flux  variation 
generally  occurs  near  L  ~  6.  (b)  Flux  level  at  L  =  6  as  a  function  of  day. 
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Predicted  Signal  Amplitude  Variation 
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Figure  3.9:  The  predicted  NLK  signal  amplitude  variation  as  observed  at  Fort  Yukon. 


Comparison  of  VLF  amplitude  data  with  GOES-7  and  SAMPEX  data  on  rel¬ 
ativistic  electron  flux  levels  show  a  clear  association  between  the  two  data  sets  in 
the  case  of  two  successive  relativistic  electron  enhancement  episodes  for  which  VLF 
data  is  available.  The  VLF  amplitude  and  particle  flux  levels  measured  on  GOES-7 
show  the  same  27  day  cycle  and  2-3  days  of  rise  and  fall  times  for  a  characteristic 
relativistic  precipitation  enhancement  event.  VLF  signal  amplitudes  exhibited  >  9 
dB  decreases  associated  with  the  electron  flux  level  enhancements  indicating  that  the 
nighttime  electron  density  at  40-70  km  altitudes  is  strongly  influenced  by  the  solar 
rotation,  via  the  relativistic  electron  enhancement  events  driven  by  the  soleir  wind. 
The  ionospheric  effect  of  the  relativistic  electron  enhancements  was  observed  only 
when  the  flux  was  above  3  x  10^  el/cm^-sr-s,  apparently  because  of  the  fact  that  the 
VLF  signature  of  the  enhancement  for  lower  fluxes  is  suppressed  by  other  ionospheric 
variations. 

Comparison  of  our  VLF  observations  with  theoretical  predictions  of  amplitude 
decreases  of  >  7  dB  obtained  using  propagation  model  calculations  provides  satisfac¬ 
tory  agreement.  Calculations  also  show  that  the  amplitude  change  associated  with  the 
lower  peaks  of  the  relativistic  electron  enhancement  are  less  than  1  dB,  not  observable 
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in  the  presence  of  larger  ionospheric  variations  associated  with  auroral  effects. 

We  conclude  on  the  basis  of  both  observation  and  the  theoretical  analysis  pre¬ 
sented  here  that  the  conductivity  of  the  nighttime  lower  ionosphere  at  subauroral 
latitudes  is  strongly  modulated  by  the  relativistic  electron  precipitation  which  accom¬ 
panies  relativistic  electron  enhancements.  High-energy  precipitation  causes  electron 
density  enhancements  in  the  D  region  of  the  ionosphere,  which  in  turn  affect  VLF 
waves  propagating  in  the  perturbed  earth-ionosphere  waveguide.  This  redization  also 
provides  the  first  evidence  of  a  detectable  influence  on  the  nighttime  lower  ionosphere 
of  solar  rotation,  imposing  a  27-day  cycle  on  top  of  other  variations  of  this  region  of 
our  atmosphere. 

Our  results  further  indicate  that  VLF  remote  sensing  can  be  a  powerful  tool 
for  investigation  of  relativistic  electron  flux  enhancements  and  their  ionospheric  and 
mesospheric  effects.  A  system  with  multiple  receiving  stations  observing  VLF  signals 
that  cross  the  affected  regions  (see  Figure  3.4)  could  be  used  to  assess  the  spatial 
distribution  of  precipitation  as  well  as  the  precipitation  flux  levels. 


3.2  Determination  of  the  Electron  Density  Profile 

The  close  association  between  relativistic  electron  flux  and  VLF  signal  level  as  re¬ 
ported  by  Demirkol  et  al.  [1999]  and  described  in  the  previous  section,  and  the 
relatively  well  defined  latitude  range  over  which  these  enhancements  occur,  form  the 
basis  of  our  method  for  determination  of  the  electron  density  profile  in  different  re¬ 
gions  along  VLF  subauroral  propagation  path. 

Noting  that  the  largest  day  to  day  variations  occur  at  Xci6,  we  model  the  geo¬ 
magnetic  latitude  distributions  of  the  precipitating  electron  fluxes  as  a  Gaussian  with 
a  center  at  Z=6  as  illustrated,  in  Figure  3.10.  The  Gaussian  spatial  distribution  of  the 
flux  level  is  illustrated  in  Figure  3.10-a.  A  stepwise  approximation  to  the  Gaussian  is 
made  as  illustrated.  Each  level  of  flux  denoted  by  a,  b,  c  and  d  results  in  an  associated 
electron  density  profile  as  shown  in  Figure  3.10-b.  The  peak  value  and  the  horizontal 
width  of  this  Gaussian  function  is  yet  to  be  determined.  In  order  to  find  the  width 
and  the  peak  of  the  Gaussian,  and  thus  to  be  able  to  identify  the  electron  density 


79 


CHAPTER  3.  AMBIENT  ELECTRON  DENSITY  PROFILE 


60 


Figure  3.10;  The  flux  level  is  modeled  as  a  Gaussian  centered  at  L  =  6.  Since  each  flux 
level  corresponds  to  a  different  electron  density  profile,  it  is  possible  to  segment  the 
propagation  path,  each  segment  including  the  corresponding  electron  density  profile 
as  shown  a,  b,  c  and  d  above. 
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Figure  3.11:  The  segmentation  of  the  VLF  propagation  path  from  NAA  to  Wasilla 
according  to  the  Gaussian  spatial  flux  model  centered  at  L  =  6. 
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profiles  corresponding  to  each  of  the  flux  levels,  the  measured  VLF  signal  amplitudes 
are  used.  For  a  set  of  different  peaks  and  widths,  a  segmentation  of  the  propagation 
path  is  made,  determined  by  the  width  of  the  Gaussian.  Figure  3.11  shows  the 
segmentation  of  the  VLF  propagation  path  from  NAA  to  Wasilla,  a  receiver  site  used 
in  the  March  1999  HAARP  Heating  campaign,  according  to  this  Gaussian  spatial  dis¬ 
tribution  model  for  the  relativistic  electron  flux.  In  each  segment  of  the  propagation 
path,  the  flux  level  is  different.  As  the  VLF  signal  propagates  from  the  transmitter  to 
the  receiver,  it  enters  the  region  affected  by  the  relativistic  electron  precipitation.  As 
the  signal  propagates  toward  higher  latitudes,  it  encounters  increasing  flux  levels,  as 
in  each  subsequent  segment  the  electron  density  in  the  ionosphere  increases  further 
and  further,  reaching  its  maximum  near  to  L~6.  At  the  center  of  the  Gaussian  distri¬ 
bution  corresponding  to  L=6,  the  flux  level  has  the  highest  value.  In  the  2D  modeling 
of  the  propagation  of  the  VLF  signal,  the  electron  density  profile  in  the  ionosphere 
over  the  region  drawn  in  blue  in  Figure  3.11  is  the  profile  produced  by  the  highest 
flux  as  shown  in  Figure  3.10.  As  the  VLF  signal  continues  to  propagate  (towards 
lower  latitudes)  the  flux  level  decreases  and  thus  the  electron  density  profile  in  the 
ionosphere  changes.  This  latitude  dependence  of  the  electron  density  is  accounted 
for  in  LWPC  by  using  different  electron  density  profiles  in  different  segments  of  the 
propagation  path. 

In  this  manner  LWPC  models  are  generated  for  a  range  of  Gaussian  models  with 
different  width  and  peak  values,  with  the  amplitudes  of  the  VLF  signal  at  the  receiver 
corresponding  to  each  of  the  different  models  calculated  as  a  byproduct  of  the  LWPC 
calculations.  Thus  it  is  possible  to  construct  a  table  where  each  Gaussian  width  and 
peak  flux  level  corresponds  to  a  particular  VLF  amplitude.  Having  such  a  table  then 
enables  us  to  identify  a  unique  pair  of  width  and  peak  by  matching  VLF  amplitudes 
at  three  frequencies,  necessarily  arriving  at  the  receiver  over  three  different  paths 
with  the  same  flux  model  imposed  on  them.  This  method  is  used  to  determine  the 
ambient  electron  density  profile  along  the  propagation  path  of  the  VLF  signal  prior 
to  the  inversion  algorithm. 
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Chapter  4 

Modified  Collision  Frequency 


Powerful  ionospheric  heaters  are  typically  designed  to  focus  an  intense  beam  of  HF 
radiation  onto  a  relatively  small  spot  with  transverse  dimensions  of  ~10  to  20  km 
at  D  region  altitudes  of  interest  here.  The  primary  physical  effect  of  a  powerful 
HF  electromagnetic  wave  in  the  D  region  is  to  heat  the  electrons  via  absorption 
of  the  wave,  which  in  turn  enhances  the  frequency  of  the  collisions  between  these 
electrons  and  the  ambient  neutrals  and  heavy  ions.  In  addition,  the  enhanced  electron 
temperature  changes  the  rates  of  chemical  reactions  which  determine  the  electron 
density  profile.  Both  electron  density  and  collision  frequency  changes  result  in  the 
modification  of  the  electrical  conductivity  of  the  heated  ionospheric  region.  Such 
a  patch  of  ionosphere  with  modified  collisional  frequency  and  electron  density  in 
turn  diffract  VLF  radio  waves  propagating  in  its  vicinity  [Barr  et  at,  1985].  Such 
diffraction  of  the  VLF  radio  waves  due  to  ionospheric  heating  by  the  HIPAS  facility 
was  illustrated  in  Section  1.4. 

Of  the  two  effect  discussed  above,  the  electron  density  modulation  resulting  from 
the  temperature  dependence  of  the  dissociative  recombination  rates  of  Oj  and  NO"*" 
[Rodriguez  and  Inan,  1994]  is  not  significant  in  the  context  of  this  work.  That  this 
is  the  case  can  be  seen  from  the  characteristic  time  of  electron  density  modulation 
(resulting  from  the  modulation  of  the  recombination  process)  which  is  given  by  [Stubbe 
et  at,  1982]  : 
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^  _  1  2.5  X  10^2 

”  aeffiVe  Ne{m-^)  ® 

For  a  typical  D  region,  Ne  ci  10®  m-^,  while  r„  is  of  the  order  of  thousands  of 
seconds.  Since  we  measure  the  effects  of  heating  over  time  scales  of  seconds,  emy 
change  in  electron  density  due  to  changes  in  the  recombination  rates  can  therefore  be 
neglected. 

Having  described  a  method  in  Chapter  3  for  the  determination  of  the  ambient 
electron  density  profile,  and  assuming  this  ambient  profile  to  be  constant  during  the 

heating  process,  we  now  discuss  the  HF  heater-induced  modification  of  the  collision 
firequency. 


4.1  HF  Heater-Induced  Modification  of  Electron 
Collision  Frequency 

For  heater  wave  power  densities  at  30  km  altitude  of  <  10“^  W/m^,  it  is  possible 
to  use  a  model  for  which  v  ~  where  is  the  thermal  velocity  [Inan  et  ai,  1991]. 
Under  these  conditions  the  electron  distribution  function  can  be  assumed  to  remain 
Maxwellian,  i.e.. 


where 


/(c)  =  N, 


t’th  =  l^thl  is  the  random  thermal  velocity  of  the  electrons, 

Ne,  is  the  number  density  of  the  electrons, 

me,  is  the  mass  of  an  electron, 

ks,  is  Boltzmann’s  constant,  and 

Tn,  is  the  ambient  temperature  of  the  electrons. 


(4.2) 
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In  this  case  the  modified  electron  collision  frequency  i/  can  be  written  as  [Maslin, 
1974]: 

where 


uq  is  the  ambient  electron  collision  frequency, 
r  is  the  absorbed  wave  power,  and 

G  is  the  fraction  of  electron  energy  lost  per  collision,  typically  taken  to  be  in  the 
range  of  1.3  x  10~^  to  6  x  10“^  [Hxtxley  and  Zaazou,  1949;  Fejer,  1970]. 


The  absorbed  wave  power  r(2)  at  each  altitude  can  be  calculated  by  propagatiing 
the  wave  incrementally  in  altitude.  The  absorption  of  the  two  propagating  waves 
corresponding  to  the  ordinary  and  extra-ordinary  modes  produce  the  total  heating. 
The  refractive  index  n  at  each  altitude  can  be  calculated  using  equation  (4.4): 


n^  =  l 


X{U  -  X) 


U{U  -X)-^Y^  + Sr 


where 


Sr  =  ^1^4  sin4  ^  +  Y^{U  -  Xy  cos2  tp, 


y  =  |F|  = 


cBq 


Bo  is  the  constant  magnetic  field  of  the  Earth, 


(4.4) 


e  is  the  charge  of  an  electron. 


rUe  is  the  mass  of  an  electron, 

LU  is  the  angular  frequency  of  the  heating  wave  in  radians/sec, 


U  =  l-iZ, 
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Z  =  i//u, 

X  =  Nee^/{eomu)^),  and 

Ip  is  the  angle  between  the  wave  vector  k  and  Bq. 

In  this  case,  the  absorbed  power  can  be  calculated  at  each  incremental  altitude  as 
{Pd/2)  exp  (2xor^A/i)  and  =  (Pd/2)  exp  {2xex^Az)  where  Xor  =  Ini(nor) , 
Xex  =  Im  (nex).  The  absorbed  power  in  each  altitude  increment  is  determined  sepa¬ 
rately  for  each  of  the  two  elliptically  polarized  components  (i.e.,  o-mode  and  x-mode) 
that  constitute  the  input  HF  wave,  with  the  total  heating  being  a  superposition  of 
that  due  to  each  component  [Tamnenko  et  ai,  1992], 

It  is  thus  possible  to  determine  the  modified  collision  frequency  profile,  given  the 
frequency,  power  density,  polarization  and  k  vector  orientation  of  the  pump  wave  and 
for  a  specified  ambient  electron  density  and  collision  frequency  profiles.  An  example 
effective  collision  frequency  profile  corresponding  to  heating  by  the  HIPAS  transmitter 
(Pd  =  9.3  X  10"^  Watts/m2;  also  see  Section  1.4)  is  shown  in  Figure  4.1. 

It  is  clear  from  the  above  discussion  that  the  resultant  altitude  profile  of  effective 
collision  frequency  is  a  function  of  the  ambient  electron  density  and  collision  frequency 
profiles.  The  ambient  electron  collision  frequency  is  largely  determined  by  the  density 
of  the  neutral  atmospheric  constituents  and  thus  does  not  exhibit  significant  variations 
from  night-to-night.  On  the  other  hand,  the  ambient  electron  density  in  the  nighttime 
lower  ionosphere  is  notoriously  variable  [Davies,  1990;  p.  143  and  Fig.  5.2].  This  high 
degree  of  variability  makes  theoretical  prediction  of  the  effective  collision  frequency 
profiles  highly  uncertain  and  underscores  the  need  for  diagnostic  experiments  that  can 
directly  determine  the  modified  collision  frequency  profile.  Even  if  the  ambient  density 
profile  is  known,  it  should  be  noted  that  the  HF  heating  model  described  in  this  section 
is  relative  simple  and  approximate,  for  example  in  terms  of  the  numerical  value  and 
altitude  dependence  of  the  quantity  G.  In  view  of  this,  accurate  determination  of  the 
modified  collision  frequency  profile  can  best  be  done  by  measurement  with  methods 
such  as  the  VLF  D  region  diagnostics  technique  described  in  this  dissertation. 

Based  on  the  fact  that  subionospheric  VLF  signals  have  been  shown  to  exhibit 
detectable  amplitude  and  phase  changes  in  response  to  heating  of  the  D  region  by  HF 
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Figure  4.1;  Calculated  disturbed  collision  frequency  profile  for  a  given  ambient  col¬ 
lision  frequency  and  electron  density  profile.  The  assumed  angle  between  the  wave 
normal  and  B  is  13.5®.  The  frequency  of  the  heater  wave  is  2.86  MHz,  and  the  two 
mode  combined  power  of  the  HF  transmitter  is  9.3  x  10“'*  watts/m^. 
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[Barr  et  al.,  1985;  Bell  et  ai,  1995]  or  VLF  waves  [Inan,  1990;  Inan  et  ai,  1992],  such 
measurements  can  potentially  be  used  to  infer  the  altitude  profile  of  effective  electron 
collision  frequency  in  a  modified  ionosphere.  The  development  and  application  of  an 
inversion  method  to  this  end  is  the  main  subject  of  this  dissertation.  In  the  next 
section,  we  describe  this  new  procedure  and  apply  it  to  two  hypothetical  examples. 
Application  of  the  method  to  actual  experimental  data  is  the  subject  of  Chapter  5. 


4.2  Inversion  Algorithm  to  Determine  u^z) 

The  amplitude  and  phase  changes,  AA  and  A<f>,  measured  by  VLF  receivers  in  re¬ 
sponse  to  HF  heating  sensitively  depend  on  ionospheric  parameters  as  described  in 
Chapter  2.  Having  determined  the  ambient  electron  density  profile  as  described  in 
Chapter  3,  the  only  unknown  parameter  affecting  AA  and  A<f)  is  the  disturbed  colli¬ 
sion  frequency  profile  i/(z)  in  the  heated  region. 

We  now  describe  the  inversion  technique  we  use  to  determine  u{z)  from  measured 
values  of  AA  and  A<p.  For  this  purpose,  we  use  three  VLF  signals  with  different 
frequencies  and  propagation  paths  so  as  to  have  six  measurements  at  our  disposal 
to  specify  six  degrees  of  freedom  for  the  collision  frequency  profile.  In  the  present 
work  these  six  degrees  of  freedom  are  chosen  to  be  the  values  of  the  electron  coUision 
frequency  at  six  different  altitudes  as  shown  in  Figure  4.2. 

Receiver  locations  are  chosen  such  that  the  three  VLF  signals  pass  under  the 
heated  region  as  illustrated  in  Figure  4.3.  The  fact  that  the  three  signals  have  dif¬ 
ferent  frequencies  ensures  that  the  resultant  amplitude  and  phase  changes  provide 
nonredundant  information. 

In  general,  we  can  define  a  calculated  change  vector  A  in  the  form  : 

A  =  {AA\,  A4>\,AAl,  A<f>l,  AA|,  A0|)  (4.5) 

where  the  vector  elements  are  calculated  using  the  3D  LWPC  code  with  a  given 
collision  frequency  profile.  Thus  : 

I'i,  i^e),  i  =  1, 3, 5, 
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Figure  4.2:  Collision  frequency  profile  is  modeled  with  six  different  values  of  collision 
frequency  at  six  different  altitudes  in  the  D-region.  The  altitude  range  sampled  is 
chosen  on  the  basis  of  the  sensitivity  of  the  subionospheric  VLF  signal  to  changes  at 
different  altitudes. 
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Figure  4.3:  Three  VLF  signals  with  different  frequencies  propagating  under  a  heated 
region  will  be  used  in  the  inversion  technique.  The  propagation  paths  of  the  signals 
from  the  transmitters  (T)  to  the  receivers  (R)  are  also  different,  since  in  general  each 
VLF  transmitter  operates  at  a  pre-allocated  frequency. 


Mi  =  1^2, 1^3, 1^4, 1^5,  i  =  2, 4, 6, 

where  are  the  6  parameters  that  define  the  collision  frequency 

profile  (Figure  4.2). 

Note  that  each  element  of  A  is  a  function  of  the  collision  frequency  profile  as  well 
as  other  parameters,  the  values  of  which  are  known. 

Given  the  measured  change  vector  =  (A>li,  A^i,  A>l2, A02,  Ayla,  A^a), 

where  Aj4i,  A/la,  A^la  are  the  three  measured  amplitude  changes,  and  A0i,  A<;62,  A^a 
are  the  corresponding  measured  phase  changes,  a  simple  algorithm  for  iteratively 
solving  for  the  collision  frequency  profile  u  =  (i/i,  1^2,  z/a,  1^4, 1/5,  i/g)  can  be  described  as 
follows  : 

1)  Start  with  an  initial  estimate  vector  =  (i/j ,  z/® ,  1/3,  i/°,  1/5),  where  the  super¬ 

script  0  corresponds  to  the  first  iteration  number. 

2)  For  the  m^^  iteration,  form  12  distinct  collision  frequency  profiles  based  on  the 

reference  points  : 
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±  5if  (4.6) 

where 

5up  =  aj{0,..,  ^  ,..,0)  and  j  =  l,2,..,6, 

position 

Qj  is  a  constant  which  will  be  referred  as  step  size,  and 
m  is  the  iteration  number. 

3)  For  each  of  the  12  profiles,  find  A. 

4)  For  the  new  profile  choose  +  Suj^,  k  being  the  index  which  minimizes 

the  error  function  defined  as: 


\ 


Ai  I P  +  |+ 


where  the  subscript  k  denotes  the  calculated  values  of  Aj4  and  A0,  —  (A -{■  A  A)  I  A, 

and  Ai,2,3  are  signal  quality  parameters  determined  by  the  signal-to-noise  ratio  at  each 
receiver  site.  Note  that  to  first  order  in  A.A  and  A0,  the  right-hand-side  of  (4.7)  is 
simply  the  common  mean  square  error  function.  The  particular  form  used  here  is 
chosen  for  convenience,  since  ^4^  is  one  of  the  output  parameters  of  the  3D  VLF 
propagation  code. 


5)  Repeat  steps  2  through  3  with  until  the  error  is  reduced  below  a  desired 
level,  i.e.,  until  |Afc  -  <  Sdesired- 


Figure  4.4  shows  the  flow  graph  for  this  algorithm. 

For  a  single  VLF  signal,  each  iteration  corresponds  to  finding  the  amplitude  and 
phase  change  for  12  different  collision  frequency  profiles  generated  by  ±  (5i^  for 
j  =  1,2,  ..,6  and  choosing  the  profile  whose  predicted  amplitude  and  phase  change 
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Figure  4.4:  Flow  graph  for  the  inversion  algorithm. 


minimizes  the  right-hand-side  of  (4.7),  as  illustrated  in  Figure  4.5  In  this  Figure 
the  amplitude  and  phase  changes  corresponding  to  the  collision  frequency  profiles 
generated  by  ±  for  j  =  1, 2, ..,  6  are  shown  as  stars.  The  amplitude  and  phase 
corresponding  to  the  initial  collision  frequency  profile  is  shown  as  an  open  circle. 
The  measured  amplitude  and  phase  change  (in  this  case  corresponding  to  =  0.14 
dB  and  A0  =  0.95°)  is  shown  in  the  upper  right  corner  and  is  the  desired  model 
output.  Among  the  different  12  profiles,  we  choose  the  profile  for  which  the  predicted 
amplitude  and  phase  changes  minimize  the  error  function.  The  amplitude  and  phase 
change  corresponding  to  this  model  is  circled.  In  the  next  iteration  step,  the  collision 
frequency  profile  corresponding  to  this  point  is  used  as  the  initial  profile. 

We  now  test  the  inversion  algorithm  by  using  the  3D  VLF  propagation  and  scat¬ 
tering  model  to  calculate  amplitude  and  phase  changes  for  each  of  the  three  signals 
(with  frequencies  /i,2,3  and  receiver  sites  Rixs)  for  a  prespecified  i/{z)  and  use  these 
as  the  measured  change  vector  A.  We  test  the  inversion  algorithm  with  two  different 
example  collision  frequency  profiles  as  shown  in  Figure  4.6.  These  collision  frequency 
profiles  are  used  to  generate  the  measured  change  vector  A  consisting  of  the  resulting 
amplitude  and  phase  change  predictions  (A^,  A4>)  as  shown  in  Table  4.1. 

In  order  to  start  the  inversion,  the  initial  collision  frequency  profile  is  required 
and  is  determined  b}^  selection  from  among  a  set  of  generic  collision  frequency  profiles 
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Figure  4.5:  For  a  single  VLF  signal  each  iteration  corresponds  to  finding  the  amplitude 
and  phase  change  for  12  different  collision  frequency  profiles  generated  by  £7"  ±  Sjv^ 
ioT  j  =  1, 2, ..,  6  and  choosing  the  profile  which  results  in  the  smallest  error  value. 
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VLF  signal 

[dB] 

Ad)  [deg] 

IdBI 

A<?!)^‘^[deg] 

Example  1 

Ri/fi 

0.15 

-0.26 

0.14 

-0.27 

Rilh 

-0.01 

0.10 

Rzth 

0.02 

■iliVJM 

lEEBU 

Example  2 

Rilfi 

0.24 

Rzlh 

-0.04 

0.65 

-0.04 

0.64 

^  Rzifs 

0.09 

0.89 

6.10 

0.90 

Table  4.1:  Amplitude  and  phase  changes  for  the  original  and  regenerated  collision 
frequency  profiles. 


(Figure  4.7),  generally  consistent  with  the  range  of  those  expected  on  the  basis  of 
the  simplified  theory  of  HF  heating  described  in  Section  4.1.  The  selection  of  is 
made  by  using  the  3D  VLF  code  with  the  generic  profiles  of  Figure  4.7  to  evaluate  the 
corresponding  amplitude  and  phase  changes  in  the  VLF  signal  closest  (in  terms  of 
the  error  function)  to  the  measured  values.  The  initial  profiles  so  chosen  are  shown  in 
Figure  4.6  and  give  the  closest  amplitude  and  phase  changes  to  the  measured  values 
in  both  examples. 

With  the  initial  profile  determined  from  among  a  finite  set  which  gives  the  best 
initial  agreement,  we  proceed  with  steps  2  through  4  of  the  inversion  algorithm.  The 
final  collision  frequency  profiles  obtained  after  m=20  iterations  are  shown  in  Figure 
4.6.  We  note  that  these  profiles  are  almost  identical  to  the  original  profiles,  and  that 
the  amplitude  and  phase  changes  corresponding  to  the  regenerated  profiles  (A|„  A((>‘^°) 
are  close  to  the  original  input  values  within  the  measurement  limitations  as  shown  in 
Table  4.1. 
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Figure  4.6:  Tests  of  the  inversion  algorithm,  (a)  Example  1;  (b)  Example  2 


Chapter  5 

March  1999  HAARP  Campaign 


As  mentioned  previously,  the  VLF  D  region  diagnostic  technique  described  in  Chapter 
4  was  developed  for  the  purpose  of  determining  the  collision  frequency  profile  in  a 
localized  ionospheric  patch,  as  part  of  the  overall  suite  of  diagnostic  instruments  for 
the  HF  Active  Auroral  Research  Program  (HAARP).  The  first  opportunity  to  field 
test  this  new  diagnostic  technique  became  available  in  March  1999,  where  the  partially 
completed  HAARP  heater  was  utilized  in  a  month  long  research  campaign.  In  this 
Chapter,  we  describe  this  campaign  and  in  particular  discuss  the  results  obtained 
with  the  new  VLF  D  region  diagnostic  method. 

The  HAARP/HIPAS  experimental  research  campaign,  conducted  in  Alaska  during 
March  1999,  was  a  complex  undertaking,  employing  two  powerful  HF  radio  transmit¬ 
ters  and  a  variety  of  radio,  radar,  and  optical  diagnostic  instruments,  located  at  a 
number  of  sites,  some  in  remote  areas  of  the  state.  A  large  number  and  wide  variety 
of  experiments  were  conducted  as  part  of  the  campaign,  carried  out  during  the  period 
of  9-29  March  1999.  The  HAARP  HF  transmitter  operated  on  19  separate  days,  with 
a  total  operation  time  of  approximately  270  hours. 


5.1  HAARP  Facility  Description 

Both  the  HAARP  and  HIPAS  facilities  include  powerful  HF  transmitters  and  a  vari¬ 
ety  of  radio  and  optical  diagnostic  instruments  used  to  conduct  active  experimental 
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investigations  to  characterize  the  physical  processes  that  can  be  initiated  and  con¬ 
trolled  in  local  regions  of  the  atmosphere,  ionosphere  and  space  \'ia  interactions  ■with 
high  power  radio  waves  [Kossey  and  Kennedy,  1999' 

The  present  HAARP  HF  transmitting  system,  located  in  Gakona.  Alaska,  includes 
a  phased-array  antenna,  consisting  of  48  elements,  with  crossed  dipole  antennas  driven 
individually  by  10  kW  transmitters,  resulting  in  a  maximum  total  radiated  power  of 
960  k\\ .  The  HAARP  is  yet  only  partially  complete.  WTien  fully  com.pleted  in  a  few 
vears.  the  HAARP  transmitter  will  consist  of  180  antenna  elements  with  a  maximum 
total  radiated  power  of  3.6  MW.  Figure  5.1  shows  a  small  section  of  the  HAARP  48 
element  HF  antenna  array. 


The  H.4ARP  HF  transmitter  can  be  operated  at  any  selected  and  authorized 
frequency  between  2.8  MHz  and  10  MHz.  It  is  possible  to  maximize  the  deposition 
of  the  radiated  energy  selectively  at  altitude.^  between  70  km  and  90  km  [D  region 
of  the  iono.sphere)  when  operating  at  2.8  MHz.  At  higher  frequencies  more  of  the 
energy  is  deposited  at  higher  altitudes.  .Although  H.4ARP  operates  between  2.8  MHz 
and  10  MHz,  the  ionospheric  effects  it  can  produce  span  16  decades  in  frequency, 
fronn  the  generation  of  ULF/ELF  waves  to  the  production  of  optical  emissions  in  the 
ionosphere.  Figure  5.2  lists  the.se  physical  effects  'J<o$sey  and  Kennedy.  1999'. 
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Figure  5.2:  Effects  spanning  16  decades  in  firequenc}^  produced  by  HAARP. 


The  Stanford  University  VLF  D  region  diagnostic  system  constitutes  only  one 
of  the  diagnostic  instruments  built  for  use  \\dth  the  HAARP  facility.  The  suite  of 
on-site  and  off-site  diagnostic  instruments  available  to  the  facility  play  an  important 
role  in  monitoring  the  local  geophysical  and  electromagnetic  background.  They  serve 
an  essential  role  during  active  experiments  employing  the  HF  transmitter,  providing 
knowledge  of  conditions  in  the  ionosphere  and  in  space  prior  to,  during,  and  after 
its  operation.  Data  collected  from  most  of  the  instruments  can  be  processed  and 
displayed  at  the  site,  and  can  be  provided  in  real-time  to  other  researchers  via  the 
Internet.  Table  5.1  shows  HAARP  Diagnostic  Instruments  [Kossey  and  Kennedy, 
19991. 

The  HIP.4S  facility,  located  near  Fairbanks,  Alaska,  has  been  in  operation  for  more 
than  a  decade.  Its  HF  transmitting  system  consists  of  8  crossed-dipole  antennas  on 
a  circular  field  with  a  diameter  of  208  meters.  The  array  has  a  total  radiated  power 
of  1  MW  at  2.85  MHz  or  4.53  MHz. 
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VLF  D  region  Diagnostic 

HF  vertical  incidence  sounder 

HF  SEE  receivers 

ELF/VTF/LF  receivers 

Magnetometer 

V  HF  (30  MHz)  Riometer 

Imaging  Riometer 

Spectrum  Monitor 

IR  Photometer 

HF  (2S  MHz)  radar 

V'HF  (50  MHz)  radar 

UHF  incoherent  scatter  radar 

Optical  Imager 

Rayleigh  LID.AR 

UHF  scintillation  monitors 

Table  5.1;  List  of  HAARP  Diagnostic  tools. 


5.2  VLF  Remote  System 


During  the  March  1999  campaign,  a  VLF  D  region  diagnostic  system  was  operated 
by  Stanford  University  to  determine  electron  temperature  changes  in  the  D  region 
during  active  HF  heating  experiments  at  the  HAARP  facility. 

The  V  LF  sensor  for  each  receiving  system  is  a  square  loop  antenna  with  an  area 
of  1./^  m~.  The  signal  received  by  the  antenna  is  amplified  by  a  preamplifier  and 
transmitted  down  a  cable  to  a  VTLF  line  receiver.  The  output  signal  from  the  line 
receiver  is  filtered  and  read  by  a  computer  via  an  A/D  card.  The  amplitudes  and 
phases  of  the  VLF  signals  from  the  three  transmitters  are  then  extracted  and  saved 
Johnson  et  al,  1999; .  Synchronization  in  time  with  the  other  sites  is  achieved  through 
the  use  of  a  GPS  receiver  to  time  stamp  the  sampled  data. 


The  amplitude  and  phase  of  three  subionospheric  V'LF  signals  propagating  through 
the  ionosphere  lying  over  the  HAARP  facility  were  measured  at  three  locations  in 


Alaska,  name.y  at  Healy,  Wasilla  and  Delta  Junction  as  shown  in  Figure  5.3.  V'LF 
signals  from  US  .Xavy  transmitters  in  Maine.  (X.AA),  Hawaii  (XPM);  and  Washington 
(XLK)  were  monitored.  At  each  location  the  diagnostic  system  was  sited  at  a  local 
high  school.  Table  5.2  shows  the  latitude  and  longitude  of  the  receivers.  Table  5.3 


shows  the  latitudes  and  longitudes  of  the  transmitters  toge 
and  the  power  of  the  radiated  VLF  signal. 


ther  with  the  frequency 


Figure  5.3  shows  the  location  of  the  receivers  and  transmitters  together  with  the 
Gicat  Ciiclc  Paths  Oi  eacn  V  LF  signal.  Figure  5.4  shows  the  heated  region  over 
the  receiver  sites  in  e-xpanded  deta.iI.  illustrating  that  the  locations  of 
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Receiver  Name 

Latitude 

Longitude 

High  School  Name 

Healy 

63°  52'  22.  T'  N 

149°  00'  49.6"  W 

Healy  High  School 

Wasilla 

61°  35'  16.9"  N 

149°  25'  58.0"  W 

Wasilla  High  School 

Delta  Junction 

63°  59'  59.6"  N 

144°  44'  00.1"  W 

Healy  Lake  High  School 

Table  5.2:  List  of  VLF  receivers  used  in  the  VLF  Diagnostic  System  during  HAARP 
1999  Campaign 


Transmitter  Name 

Latitude 

Longitude 

Frequency 

(kHz) 

Power 

(kW) 

NLK  (Jim  Creek,  Washington) 

48°  12'  N 

121°  55'  W 

24.8 

850 

NAA  (Cutler,  Maine) 

44°  39'  N 

67°  17'  N 

24.0 

1000 

NPM  (Lualualei,  Hawaii) 

21°  25'  N 

158°  09'  N 

23.4 

300 

Table  5.3:  List  of  VLF  transmitters  used  during  HAARP  1999  Campaign 


the  receivers  are  chosen  so  that  the  propagation  path  of  the  VLF  signals  from  their 
respective  transmitters  pass  within  50  km  of  the  heated  region  over  HAARP. 


5.3  Experimental  Results 

As  part  of  the  March  1999  HAARP  campaign,  15  minute  nightly  periods  were  allo¬ 
cated  for  the  VLF  D  region  diagnostic  experiment.  For  this  purpose,  the  HAARP 
HF  output  transmissions  were  100%  amplitude  modulated  with  25  Hz  square  waves 
(i.e.,  20  ms  ON/20  ms  OFF).  Figure  5.5  shows  an  example  of  the  resulting  amplitude 
changes  produced  in  the  signals  from  NPM,  NAA  and  NLK  during  the  HF  heating 
period.  The  data  shown  has  been  averaged  through  a  superposed  epoch  technique. 
The  right  hand  panels  show  the  25  Hz  modulation  of  VLF  signals  from  NPM,  NAA 
and  NLK  in  the  time  domain.  The  left  hand  panels  show  the  peaks  at  25  Hz  in 
the  Fourier  spectrum  of  the  high  resolution  data.  Figure  5.6  shows  an  example  of 
amplitude  and  phase  changes  in  the  signal  from  NLK  observed  at  Healy.  The  25  Hz 
modulation  is  clearly  apparent  both  in  the  time  domain  and  in  the  frequency  spectra. 
The  peak  at  39  Hz  in  the  NPM-Delta  Junction  spectra  is  due  to  local  interference, 
present  all  the  time,  regardless  of  whether  HAARP  HF  transmitter  is  ON. 

The  direct  association  of  the  observed  25  Hz  modulation  of  the  VLF  amplitude 
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Figure  5.3:  VLF  transmitter  and  receiver  locations. 
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Figure  5.4:  Receiver  locations  with  respect  to  the  circular  region  heated  by  HAARR 
The  receivers  are  located  such  that  the  propagation  paths  of  the  VLF  signals  pass 
through  the  heated  region. 


and  phase  with  HAARP  transmissions  is  further  illustrated  in  Figure  5.7  where  the 
spectrum  of  the  VLF  amplitude  is  plotted  for  5  sequential  5-min  periods,  the  first 
being  before  the  25  Hz  modulation  of  HAARP  is  turned  ON,  the  next  three  being 
during  which  HAARP  is  ON  and  the  last  being  immediately  after  HAARP  is  turned 
OFF.  It  is  clear  that  25  Hz  is  only  present  when  HAARP  is  ON  and  that  the  level  of 
transferred  modulation  remains  approximately  constant  during  the  15  minute  period. 
The  HAARP  VLF  D  region  diagnostic  format  was  transmitted  nearly  every  night 
during  the  March  1999  campaign,  and  VLF  data  was  acquired  at  all  three  sites, 
almost  every  night.  The  results  of  the  VLF  measurements  are  summarized  in  Tables 
5.4  through  5.6.  In  each  Table  are  shown  the  measured  absolute  amplitude  A  (average 
value  over  15-min  period  in  dB  above  1  //V/m)  of  the  corresponding  VLF  signal, 
the  amplitude  change  AA  (in  dB  above  1  jiY/m),  the  signal-to-noise  (SNRamp)  ratio 
defined  as  the  ratio  of  the  value  of  the  spectral  density  at  25  Hz  to  the  average  value 
of  spectral  density  in  the  15-35  Hz  range,  the  phase  change  A0  (in  degrees),  and  the 
phase  signal-to-noise  (SNRphs)  defined  as  the  ratio  of  the  value  of  the  spectral  density 
of  the  phase  at  25  Hz  to  its  average  value  in  the  15-35  Hz  range.  The  Table  for  the 
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NPM  at  Delta  Junction 


NAA  at  Wasilla 


NLK  at  Healy 


Time  (s) 


FFT  of  Amplitude  data 


1 

0.5 


10  20  30  40 

Frequency  (Hz) 


Figure  5.5:  Superposed  epoch  analysis  results  and  spectra  of  the  NPM,  NAA  and 
NLK  signal  amplitudes  observed  at  Delta  Junction,  Wasilla  and  Healy  respectively' 
between  05:45-06:00  UT  on  03/25/1999.  The  25  Hz  modulation  is  apparent  on  the 
NLK  and  NPM  signals.  The  peak  at  39  Hz  on  the  NPM  signal  is  due  to  local 
interference. 
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03/15/99  05:45-06:00  UT 


0  0.2  0.4  0.6  0.8  1 

Time  (s) 


FFT  of  amnlitude  data 


Frequency  (Hz) 


Figure  5.6:  Amplitude  and  phase  modulations  of  the  VLF  signal  from  NLK  as  received 
at  Healy.  A  25  Hz  modulation  is  clear  both  in  the  time  and  frequency  domains. 
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NPM  signal  measured  at  Delta  Junction  is  short  since  the  NPM  VLF  transmitter  in 
Lualuelai.  Hawaii  was  undergoing  maintenance  (and  hence  not  operational)  through 
most  of  the  March  1999  campaign. 

The  HAARP  VLF  D  region  diagnostic  system  was  originally  designed  for  the  fully 
developed  HAARP  HF  heater,  planned  to  be  operating  at  a  total  radiated  power 
level  of  3.6  MW.  Accordingly,  the  range  of  results  obtained  (Tables  5.4-5.6)  during 
the  March  1999  campaign,  conducted  at  a  radiation  power  level  of  960  kW,  are  indeed 
quite  encouraging,  even  though  the  signal-to-noise  ratios  were  marginal  (SNR  ~  1) 
in  a  number  of  days  and  were  too  low  (SNR  <  1)  in  others.  Unfortunately,  for  many 
of  the  days  during  the  campaign  the  VLF  transmitter  NPM  was  not  operational,  as 
mentioned  above.  On  days  when  NPM  was  operating,  the  SNR  at  one  or  more  of  the 
stations  was  too  low.  Thus,  it  was  not  possible  to  measure  detectable  amplitude  and 
phase  changes  simultaneously  on  VLF  signals  from  three  different  propagation  paths. 
Accordingly,  the  range  of  data  available  only  allowed  application  testing  of  the  VLF 
diagnostic  in  a  case  involving  simultaneous  observation  of  two  identifiable  amplitude 
and  phase  changes  on  two  different  signals.  The  application  of  the  inversion  algorithm 
to  the  VLF  observations  for  this  case  is  described  in  the  next  section. 

5.4  Example  Inversion  of  VLF  data  from  25  March 
1999 

The  last  3  days  of  the  campaign  are  the  only  days  with  VLF  amplitude  and  phase 
data  from  all  of  the  3  VLF  transmitters,  but  as  mentioned  above,  the  SNR  of  the 
NAA  signal  was  too  low.  We  now  apply  the  VLF  D  region  diagnostic  method  to  one 
of  the  few  days  ,  namely  3/25/99,  for  which  amplitude/phase  changes  were  observed 
(with  SNR>  1)  on  two  VLF  signals  in  order  to  present  a  step-by-step  description  of 
the  inversion  procedure. 
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Figure  5.7:  NLK  amplitude  spectrum  from  06:40  to  07:05  in  5  minute  intervals  on 
03/10/1999.  It  is  clear  that  the  amplitude  modulation  is  due  to  HAARP. 
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UTC 

Date 

UTC 

Time 

NLK  A 

(dB  above  l/xY /m) 

(dB) 

SNRamp 

A<f) 

(deg) 

3/09/99 

0645-0700 

66.75 

0.28 

25.78 

0.08 

3/10/99 

0645-0700 

68.65 

0.18 

3/11/99 

0445-0500 

70.72 

0.03 

5.44 

■iKa 

3/15/99 

0545-0600 

70.65 

0.16 

53.00 

13.53 

3/16/99 

0645-0700 

67.01 

0.17 

13.16 

nil 

21.42 

3/17/99 

0645-0700 

69.57 

0.10 

16.81 

0.48 

11.73 

3/18/99 

0645-0700 

72.18 

0.18 

27.28 

1.64 

2.27 

3/19/99 

0545-0600 

67.95 

0.04 

2.98 

1.40 

16.45 

3/20/99 

0545-0600 

^  67.73 

0.05 

8.13 

0.32 

7.66 

3/22/99 

0545-0600 

69.11 

0.04 

6.25 

1.08 

28.66 

3/23/99 

0545-0600 

70.56 

0.08 

10.55 

0.58 

10.72 

3/24/99 

0645-0700 

No  Data 

- 

- 

- 

- 

03/25/99 

0545-0600 

66.83 

0.14 

0.96 

8.50 

03/26/99 

0745-0800 

70.19 

0.00 

0.52 

6.04 

0.90 

03/27/99 

0845-0900 

72.69 

i03 

0.40 

20.98 

Table  5.4:  List  of  the  amplitude  and  phase  measurements  at  Healy. 


UTC 

Date 

UTC 

Time 

NAA  A 

(dB  above  l//V/m) 

AA 

(dB) 

SNRamp 

A<p 

(deg) 

03/09/99 

0645-0700 

No  Data 

- 

- 

HDH 

03/1099 

0645-0700 

■fiBai 

iKiiTta 

03/11/99 

0445-0500 

67.21 

EE3 

03/15/99 

0545-0600 

66.13 

Eia 

beeh 

03/16/99 

0645-0700 

67.18 

Kini 

BiM 

03/17/99 

0645-0700 

72.52 

0.08 

03/18/99 

0645-0700 

72.41 

0.12 

BiH 

03/19/99 

0545-0600 

69.96 

1.21 

0.10 

03/20/99 

0545-0600 

62T6 

EEa 

2.14 

0.24 

BlBi 

03/22/99 

0545-0600 

72.75 

iilBai 

3.17 

0.08 

WObkm 

03/23/99 

0545-0600 

66.64 

iww 

0.56 

0.08 

03/24/99 

0645-0700 

72.58 

0.00 

0.19 

HBM 

bbh 

03/25/99 

0545-0600 

69.73 

0.54 

OEI 

03/26/99 

0745-0800 

74.65 

1.23 

0.48 

03/27/99 

0845-0900 

0.49 

1.38 

Table  5.5:  List  of  the  amplitude  and  phase  measurements  at  Wasilla. 
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UTC 

Date 

UTC 

Time 

NPM  A 

(dB  above  1/zV/m) 

AA 

(dB) 

SNRamp 

A0 

(deg) 

SNRphs 

03/09-24/99 

Various 

Transmitter  OFF 

- 

- 

03/25/99 

0545-0600 

69.78 

-7.42 

0.87 

03/26/99 

0745-0800 

71.51 

0.02 

0.75 

03/27/99 

0845-0900 

70.81 

0.06 

1.33 

Table  5.6:  List  of  the  amplitude  and  phase  measurements  at  Delta  Junction. 


5.4.1  Data  Selection  and  Calibration 

The  VLF  amplitude  and  phase  data  from  NLK  and  NPM  received  on  3/25/99  re¬ 
spectively  at  Delta  Junction  and  Healy  are  shown  in  Figure  5.8.  The  NAA  signal  at 
Wasilla  did  not  exhibit  any  detectable  changes  as  indicated  in  Table  5.5.  We  see  that 
both  NLK-Healy  and  NPM-Delta  Junction  signals  clearly  show  the  25  Hz  modula¬ 
tion.  (NLK  amplitude  shows  another  peak  at  ~  39  Hz,  which  is  due  to  local  noise  and 
is  visible  during  other  periods  when  the  HAARP  transmitter  was  not  operational.) 

For  comparison,  the  spectra  of  both  raw  data  and  the  superposed  epoch  analysis 
for  the  NAA-Wasilla  signal  are  shown  in  Figure  5.9.  There  is  no  clear  peak  at  25  Hz 
showing  any  indication  of  HAARP  modulation. 

A  block  diagram  representation  of  the  VLF  receiver  at  each  site  is  shown  in  Figure 
5.10.  In  order  to  compare  the  recorded  VLF  values  with  LWPC  calculations  (which 
give  the  value  of  the  vertical  electric  field  measured  at  the  receivers),  it  is  necessary 
to  determine  the  amplitude  of  the  vertical  electric  field  measured  at  each  site  in 
fiV/m.  In  our  experiment,  we  utilize  vertically  deployed  magnetic  loop  antennas 
which  measure  the  horizontal  magnetic  field,  simply  related  to  the  electric  field,  i.e. 
Ejn  =  cBin,  where  c  is  the  speed  of  light.  The  measured  units  generated  by  the 
receiver  (Dout)  is  eventually  expressed  in  terms  of  the  equivalent  field  strength  (jFin). 

The  overall  transfer  function,  Dont/  Em  consists  of  three  gain  stages: 

•  Vin/Ein,  the  gain  of  the  antenna  with  respect  to  the  effective  input  electric  field 
strength.  The  effective  field  Ei^  has  the  value  Ein  =  cRjn,  where  3^  is  the 
magnetic  field  measured  by  the  loop.  The  gain  of  the  air-core  loop  antenna 
used  in  the  Stanford  VLF  receivers  is  ~1.88. 
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NPM  at  Delta  Junction 


FFT  of  Amplitude  data 


10  20  30  40 


FFT  of  Phase  data 
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Figiure  5.8:  On  the  left  panels,  superposed  epoch  analysis  of  amplitude  and  phase  data 

from  NLK  and  NPM  are  sho\\m.  On  the  right  panels,  the  spectra  of  the  respective 
data  are  shown. 
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FFT  of  Superposed  Phase  data 


Frequency,  [Hz] 


Figure  5.9:  NAA  amplitude  shows  no  25  Hz  modulation  on  03/25/1999. 


Figure  5.10:  The  overall  transfer  function  Pout/  ^in  consists  of  three  transfer  func¬ 
tions,  gain  of  the  antenna,  gain  of  the  analog  front  end,  and  gain  of  the  distal  receiver. 
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•  K)ut/Mn)  the  gain  of  the  analog  front  end.  This  gain  depends  on  the  frequency  of 
the  VLF  signal  being  received.  In  order  to  determine  the  analog  front  end  gain,  a 
dummy  loop  consisting  of  a  lumped  circuit  with  terminal  impedance  identical  to 
the  input  impedance  of  the  loop  antenna  is  connected  to  the  antenna  terminals. 

A  voltage  VJn  is  then  injected  by  a  function  generator,  and  the  voltage  at  the 
output  of  the  analog  front  end  is  measured  at  selected  frequencies.  From  these 
measurements  plots  of  the  analog  gain  vs  frequency  are  generated,  which  can 
be  used  to  find  the  gain  at  any  frequency  by  interpolation.  Figure  5.11  shows 
this  calibration  plots  for  3  different  receivers. 

•  Doxit/Vout,  the  gain  of  the  digital  receiver.  This  gain  mainly  depends  on  the 
dynamic  range  of  the  A/D  card  used  and  the  representation  of  the  digital  num¬ 
bers.  In  our  case,  the  A/D  card  has  an  input  range  of  ±5V  and  the  values  are 
recorded  as  16  bit  unsigned  integers,  which  gives  a  digital  gain  of  2^®/10. 

Further  details  on  the  operation  and  calibration  of  the  narrowband  system  used 
in  the  Stanford  VLF  receivers  can  be  found  at  http://ww-star.stanford.edu/'hail/. 

5.4.2  The  Ambient  Electron  Density  Profile 

Once  the  measured  VLF  amplitudes  are  converted  to  field  values  in  //V/m,  the  next 
step  is  to  determine  the  ambient  electron  density  profile  for  which  2D  LWPC  calcu¬ 
lations  as  described  in  Chapter  3  produce  approximately  the  same  measured  values. 
For  this  step  the  collision  frequency  profile  is  assumed  to  be  a  constant  standard 
profile  (Figure  1.3-b)  throughout  the  propagation  path  for  all  3  different  receivers, 
and  the  electron  density  profile  is  assumed  to  be  described  by  a  Gaussian  flux  model 
described  in  Chapter  3.  The  determination  of  the  ambient  electron  density  profile 
relies  upon  a  stored  data  file  which  contains  the  absolute  (/xV/m)  values  calculated 
with  the  2D  LWPC  code  for  different  models  which  is  used  to  select  the  best  fit  to  the 
measured  values  as  the  ambient  electron  density  profile.  A  computer  program  written 
for  this  purpose  reads  in  the  data  file,  finds  the  best  fit,  and  automatically  generates 
the  LWPC  model  input  files  (see  Section  2.1)  for  three  different  propagation  paths 
for  the  next  steps.  Figure  5.13  illustrates  the  segmentation  of  the  propagation  path 
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Figure  5.11:  Calibration  plots  for  the  receivers  at  Wasilla,  Delta  Junction  and  Healy. 
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into  segments  having  different  electron  density  profiles.  Figure  5.12  shows  the  electron 
density  profiles  that  are  used  in  generating  the  stored  data  file. 


Figure  5.12;  Different  electron  density  profiles  used  in  the  modeling  of  the  Gaussian 
spatial  distribution  centered  at  L  =  6. 


Following  the  procedure  described  above,  and  using  the  field  values  (in  fiV/m) 
measured  at  three  different  sites,  the  best  Gaussian  model  for  the  electron  density  is 
found.  The  measured  field  values  (normalized  by  the  power  levels  of  the  transmitters) 
and  the  values  produced  by  2D  LWPC  code  for  the  electron  density  profile  determined 
in  this  was  are  given  in  Table  5.7.  Note  that  the  NAA-Wasilla  data  was  not  included 
in  this  determination  since  the  signal-to-noise  ratio  for  this  case  was  less  than  unity. 

The  electron  density  profiles  corresponding  to  the  different  segments  shown  in  ^ 

Figure  5.13  are  plotted  in  Figure  5.14.  As  indicated  in  Chapter  4,  the  electron  density 
profile  in  the  D  region  of  the  ionosphere  can  be  assumed  to  be  constant  during 
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Figure  5.13;  The  GCP  paths  of  the  three  VLF  signals  are  divided  into  segments 
where  the  electron  density  profiles  are  different  and  correspond  to  a  Gaussian  fiux 
level  model  as  described  in  Chapter  3. 


Measured  field  value 

(dB  above  IfiY/m) 

Gaussian  Model  Result 

(dB  above  1/iV/m) 

NLK 

34.27 

NPM 

31.61 

31.43 

Table  5.7:  The  measured  values  of  the  transmitter  signals  at  the  receivers  and  model 
results. 
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HAARP  modulation.  Thus  the  profiles  shown  in  Figure  5.14  are  used  in  the  later  steps 
and  specifically  in  the  next  section. 

5.4.3  The  Modified  Collision  Frequency  Profile 

The  disturbed  collision  frequency  profile  is  determined  using  the  inversion  algorithm 
explained  in  Chapter  4.  It  is  first  necessary  to  identify  an  initial  guess  for  the  six 
parameters  that  describe  the  disturbed  collision  frequency  profile.  For  this  purpose, 
we  select  the  one  profile  (from  among  a  prescribed  set  as  described  in  Chapter  4) 
which  gives  the  amplitude  and  phase  changes  closest  to  the  measured  ones.  The 
initial  profile  selected  in  this  way  is  indicated  in  Figure  5.16. 

With  the  initial  profile  so  determined,  the  final  step  is  to  use  the  inversion  algo¬ 
rithm  described  in  Chapter  4.  The  iterations  for  the  NLK  amplitude  and  phase  are 
shown  in  Figure  5.15  in  the  same  format  as  in  Figure  4.5.  The  final  collision  frequency 
profile  is  shown  in  Figure  5.16  as  a  solid  line. 

It  is  interesting  to  note  that  the  final  profile  that  we  arrive  at  with  the  inversion 
method  is  substantially  different  from  the  theoretical  profile  calculated  using  the 
simple  theory  of  HF  heating  described  in  Section  4.1.  The  inferred  profile  (i.e.,  that 
which  is  marked  ’’final”  in  Figure  5.16)  indicates  that  the  electron  collision  firequency 
is  enhanced  by  nearly  an  order  of  magnitude  at  90  km  altitude  and  by  nearly  a 
factor  of  30  at  80  km  altitude.  The  sharp  variations  in  the  profile  below  ~75  km 
altitude  are  probably  not  real,  and  may  be  an  artifact  of  the  fact  that  we  have  used 
only  four  measurements  (i.e.,  two  VLF  signals)  to  characterize  a  6-parameter  profile. 
This  observation  underscores  the  need  to  use  at  least  6-parameters  to  characterize 
the  disturbed  collision  frequency  profile,  as  the  Stanford  VLF  D  region  diagnostic 
system  is  designed  for.  Nevertheless  the  maximum  heating  does  appear  to  occur  near 
80  km  altitude,  a  result  which  may  be  a  manifestation  of  the  HAARP  HF  frequency 
used  for  heating  on  this  day  (2.87  MHz).  Follow-on  heating  experiments  at  other 
HF  frequencies  may  further  clarify  the  manner  in  which  HF  heater  energy  is  actually 
deposited  in  the  lower  ionosphere. 
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Figure  5.16:  Collision  frequency  profile  used  as  an  initial  estimation,  the  final  profile 
generated  at  the  end  of  the  inversion,  and  the  theoretically  expected  profile  is  shown. 
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Chapter  6 


Summary  and  Suggestions  for 
Future  Work 


6.1  Summary 

We  have  developed  a  new  method  for  determination  of  the  electron  collision  frequency 
profile  within  a  localized  region  of  the  ionosphere  modified  in  a  transient  manner  by 
radiowave  heating,  via  inversion  of  measurements  of  associated  subionospheric  VLF 
signal  changes.  Localized  disturbances  in  the  D  region  of  the  ionosphere  are  routinely 
produced  by  powerful  HF  transmitting  facilities  such  as  the  High  Power  Auroral  Stim¬ 
ulation  (HIPAS)  facility  and  the  High  Frequency  Active  Auroral  Research  Program 
(HAARP),  as  well  as  by  lightning  discharges.  The  new  VLF  diagnostic  method  which 
is  the  subject  of  this  dissertation  takes  advantage  of  the  fact  that  electron  density  and 
temperature  changes  in  the  D  region  are  sensitively  manifested  in  changes  in  the  am¬ 
plitude  and  phase  of  subionospheric  Very  Low  Frequency  (VLF)  signals  propagating 
beneath  the  perturbed  region. 

Amplitude  and  phase  changes  of  VLF  waveguide  signals  scattered  from  artificially 
heated  ionospheric  patches  can  be  measured  with  appropriately  placed  VLF  receivers. 
The  new  inversion  algorithm  is  designed  to  determine  the  altitude  profile  of  electron 
collision  frequency  within  the  patch  by  using  the  measured  amplitude  and  phase  of 
three  different  VLF  signals  at  three  separate  frequencies.  The  algorithm  is  primarily 
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based  on  the  recursive  application  of  the  three  dimensional  version  of  the  Long  Wave 
’•  Propagation  Capability  (LWPC)  code  used  to  model  the  subionospheric  propagation 

and  scattering  of  VLF  signals  in  the  earth-ionosphere  waveguide  in  the  presence  of 
ionospheric  disturbances.  The  amplitude  and  phase  changes  at  each  receiver  produced 
by  different  collision  frequency  profiles  were  found  to  be  detectably  different  from  each 
other  even  though  the  measured  signals  spanned  a  narrow  range  of  frequency  (24-25 
kHz) ,  and  the  inversion  algorithm  was  successfully  utilized  to  determine  the  modified 
collision  frequency  profile. 

The  new  algorithm  is  applied  to  experimental  data  acquired  during  March  1999, 
when  a  VLF  D  region  diagnostic  system  was  operated  in  Alaska  in  association  with  an 
ionospheric  heating  campaign.  The  aim  of  the  experiment  was  to  determine  electron 
collision  frequency  and  temperature  changes  in  the  D  region,  during  active  HF  heating 
experiments  conducted  with  the  HAARP  facility.  For  this  purpose,  the  HF  (2.87 
MHz)  transmissions  from  the  HAARP  system  were  100%  amplitude  modulated  with  a 
25  Hz  square  wave  (ON/OFF)  pattern,  for  15-minute  periods  each  night.  The  transfer 
of  this  modulation  to  the  VLF  signals  passing  through  the  heated  region  was  detected 
through  amplitude  and  phase  changes  produced  in  the  signal  from  NPM,  NAA  and 
NLK,  using  a  superposed  epoch  analysis  of  data  acquired  during  the  same  15  minutes. 
The  HAARP  25  Hz  modulation  was  clearly  apparent  both  in  the  time  domain  and 
frequency  domain  for  most  of  the  days  on  at  least  one  of  the  NPM,  NLK  or  NAA 
transmitter  signals.  The  inversion  algorithm  for  both  finding  the  ambient  electron 
density  profile  and  the  modified  collision  fi'equency  profile  was  applied  successfully  to 
the  available  data. 

The  altitude  profile  of  electron  collision  frequency  produced  by  an  HF  ionospheric 
heater  depends  sensitively  on  the  electron  density  profile  in  the  lower  ionosphere.  For 
this  purpose,  we  have  shown  that  it  is  possible  to  take  advantage  of  a  newly  observed 
association  between  relativistic  electron  enhancements  and  nighttime  D  region  elec¬ 
tron  density  [Demirkol  et  al,  1999],  to  model  propagation  of  VLF  waves  in  subauroral 
/  regions,  and  to  assess  the  ambient  electron  density  profile  from  measured  values  of 

VLF  signal  amplitudes.  We  were  thus  able  to  determine  the  ambient  electron  den- 
\  sity  profile  along  the  propagation  path  of  the  VLF  signals  for  the  HAARP  March 
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1999  experiment,  prior  to  applying  the  inversion  algorithm  for  determination  of  the 
collision  frequency  profile. 

6.2  Suggestions  for  Future  Work 

6.2.1  Improved  Electron  Density  Model 

Determination  of  the  ambient  electron  density  profile  is  a  necessary  step  for  inferring 
the  modified  collision  frequency  profile  in  a  disturbed  region  of  ionosphere  by  an  HF 
heater.  In  this  work,  the  ambient  electron  density  profile  on  the  propagation  path  of 
the  VLF  signal  was  assumed  to  be  modulated  by  relativistic  electron  precipitation  flux 
level.  Since  the  latitude  range  that  was  of  interest  was  in  the  subauroral  regions  where 
the  precipitation  events  frequently  take  place,  the  Gaussian  modeling  of  the  spatial 
distribution  of  the  flux  level  was  applied  successfully.  However,  in  the  determination 
of  a  Gaussian  model  which  would  best  fit  the  measured  VLF  signal  levels,  only  the 
amplitude  of  the  VLF  data  was  used  for  comparison.  Thus  2  unknown  parameters  (the 
width  and  the  peak  of  the  Gaussian)  were  determined  by  3  amplitude  measurement 
of  VLF  signal.  In  the  future,  the  additional  use  of  VLF  phase  data  may  allow  us 
to  relax  some  of  the  fixed  parameters  in  the  current  model  (such  as  the  center  of 
the  Gaussian  distribution)  allowing  for  more  flexible  modeling  of  the  electron  density 
profile.  An  alternative  method  would  be  inferring  the  ambient  electron  density  profile 
using  an  inversion  method  similar  to  described  in  Chapter  4,  with  6  measurements 
and  six  degrees  of  freedom. 

6.2.2  Future  Experiments 

Results  from  March  1999  HAARP  campaign  were  encouraging  for  the  future  utiliza¬ 
tion  of  the  D  region  VLF  diagnostic  system.  Although  the  VLF  diagnostic  system  was 
designed  for  use  with  full  power  operation  of  the  HAARP  facility,  it  was  possible  to 
observe  amplitude  and  phase  modulations  in  many  cases  with  very  reasonable  signal- 
to-noise  ratio.  When  the  HAARP  heating  facility  is  fully  functional,  measurements 
and  phase  changes  with  substantially  higher  SNR  will  be  available  from  each  of  the 


122 


CHAPTER  6.  SUMMARY  AND  SUGGESTIONS  FOR  FUTURE  WORK  103 


three  receiver  sites,  making  the  diagnostic  system  more  effective.  It  is  also  possible  to 
increase  the  number  of  received  VLF  signals  in  order  to  increase  the  number  of  free 
variables  in  the  modeling  of  the  collision  frequency  profile  and/or  to  generate  a  test 
case  for  determining  the  ambient  electron  density  profile.  For  example,  by  observing 
5  distinct  VLF  signals  which  have  propagated  over  HAARP,  we  could  allocate  5  de¬ 
grees  of  freedom  to  the  ambient  electron  density  profile  and  another  5  to  the  electron 
collision  frequency  profile.  In  this  case  it  would  not  be  necessary  to  predetermine  the 
ambient  electron  density  profile  by  other  methods,  such  as  those  described  in  Chapter 


3. 


